home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / liboctave / dMatrix.cc < prev    next >
C/C++ Source or Header  |  1997-07-27  |  64KB  |  3,295 lines
  1. // Matrix manipulations.
  2. /*
  3.  
  4. Copyright (C) 1996 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  21.  
  22. */
  23.  
  24. #if defined (__GNUG__)
  25. #pragma implementation
  26. #endif
  27.  
  28. #ifdef HAVE_CONFIG_H
  29. #include <config.h>
  30. #endif
  31.  
  32. #include <cfloat>
  33.  
  34. #include <iostream.h>
  35.  
  36. #include "byte-swap.h"
  37. #include "dbleAEPBAL.h"
  38. #include "dbleDET.h"
  39. #include "dbleSCHUR.h"
  40. #include "dbleSVD.h"
  41. #include "f77-fcn.h"
  42. #include "lo-error.h"
  43. #include "lo-ieee.h"
  44. #include "lo-mappers.h"
  45. #include "lo-utils.h"
  46. #include "mx-base.h"
  47. #include "mx-inlines.cc"
  48. #include "oct-cmplx.h"
  49.  
  50. // Fortran functions we call.
  51.  
  52. extern "C"
  53. {
  54.   int F77_FCN (dgemm, DGEMM) (const char*, const char*, const int&,
  55.                   const int&, const int&, const double&,
  56.                   const double*, const int&,
  57.                   const double*, const int&,
  58.                   const double&, double*, const int&,
  59.                   long, long);
  60.  
  61.   int F77_FCN (dgeco, DGECO) (double*, const int&, const int&, int*,
  62.                   double&, double*);
  63.  
  64.   int F77_FCN (dgesl, DGESL) (const double*, const int&, const int&,
  65.                   const int*, double*, const int&);
  66.  
  67.   int F77_FCN (dgedi, DGEDI) (double*, const int&, const int&,
  68.                   const int*, double*, double*,
  69.                   const int&);
  70.  
  71.   int F77_FCN (dgelss, DGELSS) (const int&, const int&, const int&,
  72.                 double*, const int&, double*,
  73.                 const int&, double*, double&, int&,
  74.                 double*, const int&, int&);
  75.  
  76.   // Note that the original complex fft routines were not written for
  77.   // double complex arguments.  They have been modified by adding an
  78.   // implicit double precision (a-h,o-z) statement at the beginning of
  79.   // each subroutine.
  80.  
  81.   int F77_FCN (cffti, CFFTI) (const int&, Complex*);
  82.  
  83.   int F77_FCN (cfftf, CFFTF) (const int&, Complex*, Complex*);
  84.  
  85.   int F77_FCN (cfftb, CFFTB) (const int&, Complex*, Complex*);
  86.  
  87.   int F77_FCN (dlartg, DLARTG) (const double&, const double&, double&,
  88.                 double&, double&);
  89.  
  90.   int F77_FCN (dtrsyl, DTRSYL) (const char*, const char*, const int&,
  91.                 const int&, const int&, const double*,
  92.                 const int&, const double*, const int&,
  93.                 const double*, const int&, double&,
  94.                 int&, long, long);
  95.  
  96.   double F77_FCN (dlange, DLANGE) (const char*, const int&,
  97.                    const int&, const double*,
  98.                    const int&, double*); 
  99.  
  100.   int F77_FCN (qzhes, QZHES) (const int&, const int&, double*,
  101.                   double*, const long&, double*);
  102.  
  103.   int F77_FCN (qzit, QZIT) (const int&, const int&, double*, double*,
  104.                 const double&, const long&, double*,
  105.                 int&);
  106.  
  107.   int F77_FCN (qzval, QZVAL) (const int&, const int&, double*,
  108.                   double*, double*, double*, double*,
  109.                   const long&, double*);
  110. }
  111.  
  112. // Matrix class.
  113.  
  114. Matrix::Matrix (const RowVector& rv)
  115.   : MArray2<double> (1, rv.length (), 0.0)
  116. {
  117.   for (int i = 0; i < rv.length (); i++)
  118.     elem (0, i) = rv.elem (i);
  119. }
  120.  
  121. Matrix::Matrix (const ColumnVector& cv)
  122.   : MArray2<double> (cv.length (), 1, 0.0)
  123. {
  124.   for (int i = 0; i < cv.length (); i++)
  125.     elem (i, 0) = cv.elem (i);
  126. }
  127.  
  128. Matrix::Matrix (const DiagMatrix& a)
  129.   : MArray2<double> (a.rows (), a.cols (), 0.0)
  130. {
  131.   for (int i = 0; i < a.length (); i++)
  132.     elem (i, i) = a.elem (i, i);
  133. }
  134.  
  135. // XXX FIXME XXX -- could we use a templated mixed-type copy function
  136. // here?
  137.  
  138. Matrix::Matrix (const charMatrix& a)
  139.   : MArray2<double> (a.rows (), a.cols ())
  140. {
  141.   for (int i = 0; i < a.rows (); i++)
  142.     for (int j = 0; j < a.cols (); j++)
  143.       elem (i, j) = a.elem (i, j);
  144. }
  145.  
  146. bool
  147. Matrix::operator == (const Matrix& a) const
  148. {
  149.   if (rows () != a.rows () || cols () != a.cols ())
  150.     return false;
  151.  
  152.   return equal (data (), a.data (), length ());
  153. }
  154.  
  155. bool
  156. Matrix::operator != (const Matrix& a) const
  157. {
  158.   return !(*this == a);
  159. }
  160.  
  161. Matrix&
  162. Matrix::insert (const Matrix& a, int r, int c)
  163. {
  164.   Array2<double>::insert (a, r, c);
  165.   return *this;
  166. }
  167.  
  168. Matrix&
  169. Matrix::insert (const RowVector& a, int r, int c)
  170. {
  171.   int a_len = a.length ();
  172.   if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
  173.     {
  174.       (*current_liboctave_error_handler) ("range error for insert");
  175.       return *this;
  176.     }
  177.  
  178.   for (int i = 0; i < a_len; i++)
  179.     elem (r, c+i) = a.elem (i);
  180.  
  181.   return *this;
  182. }
  183.  
  184. Matrix&
  185. Matrix::insert (const ColumnVector& a, int r, int c)
  186. {
  187.   int a_len = a.length ();
  188.   if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
  189.     {
  190.       (*current_liboctave_error_handler) ("range error for insert");
  191.       return *this;
  192.     }
  193.  
  194.   for (int i = 0; i < a_len; i++)
  195.     elem (r+i, c) = a.elem (i);
  196.  
  197.   return *this;
  198. }
  199.  
  200. Matrix&
  201. Matrix::insert (const DiagMatrix& a, int r, int c)
  202. {
  203.   int a_nr = a.rows ();
  204.   int a_nc = a.cols ();
  205.  
  206.   if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
  207.     {
  208.       (*current_liboctave_error_handler) ("range error for insert");
  209.       return *this;
  210.     }
  211.  
  212.   fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
  213.  
  214.   for (int i = 0; i < a.length (); i++)
  215.     elem (r+i, c+i) = a.elem (i, i);
  216.  
  217.   return *this;
  218. }
  219.  
  220. Matrix&
  221. Matrix::fill (double val)
  222. {
  223.   int nr = rows ();
  224.   int nc = cols ();
  225.   if (nr > 0 && nc > 0)
  226.     for (int j = 0; j < nc; j++)
  227.       for (int i = 0; i < nr; i++)
  228.     elem (i, j) = val;
  229.  
  230.   return *this;
  231. }
  232.  
  233. Matrix&
  234. Matrix::fill (double val, int r1, int c1, int r2, int c2)
  235. {
  236.   int nr = rows ();
  237.   int nc = cols ();
  238.   if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
  239.       || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
  240.     {
  241.       (*current_liboctave_error_handler) ("range error for fill");
  242.       return *this;
  243.     }
  244.  
  245.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  246.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  247.  
  248.   for (int j = c1; j <= c2; j++)
  249.     for (int i = r1; i <= r2; i++)
  250.       elem (i, j) = val;
  251.  
  252.   return *this;
  253. }
  254.  
  255. Matrix
  256. Matrix::append (const Matrix& a) const
  257. {
  258.   int nr = rows ();
  259.   int nc = cols ();
  260.   if (nr != a.rows ())
  261.     {
  262.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  263.       return Matrix ();
  264.     }
  265.  
  266.   int nc_insert = nc;
  267.   Matrix retval (nr, nc + a.cols ());
  268.   retval.insert (*this, 0, 0);
  269.   retval.insert (a, 0, nc_insert);
  270.   return retval;
  271. }
  272.  
  273. Matrix
  274. Matrix::append (const RowVector& a) const
  275. {
  276.   int nr = rows ();
  277.   int nc = cols ();
  278.   if (nr != 1)
  279.     {
  280.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  281.       return Matrix ();
  282.     }
  283.  
  284.   int nc_insert = nc;
  285.   Matrix retval (nr, nc + a.length ());
  286.   retval.insert (*this, 0, 0);
  287.   retval.insert (a, 0, nc_insert);
  288.   return retval;
  289. }
  290.  
  291. Matrix
  292. Matrix::append (const ColumnVector& a) const
  293. {
  294.   int nr = rows ();
  295.   int nc = cols ();
  296.   if (nr != a.length ())
  297.     {
  298.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  299.       return Matrix ();
  300.     }
  301.  
  302.   int nc_insert = nc;
  303.   Matrix retval (nr, nc + 1);
  304.   retval.insert (*this, 0, 0);
  305.   retval.insert (a, 0, nc_insert);
  306.   return retval;
  307. }
  308.  
  309. Matrix
  310. Matrix::append (const DiagMatrix& a) const
  311. {
  312.   int nr = rows ();
  313.   int nc = cols ();
  314.   if (nr != a.rows ())
  315.     {
  316.       (*current_liboctave_error_handler) ("row dimension mismatch for append");
  317.       return *this;
  318.     }
  319.  
  320.   int nc_insert = nc;
  321.   Matrix retval (nr, nc + a.cols ());
  322.   retval.insert (*this, 0, 0);
  323.   retval.insert (a, 0, nc_insert);
  324.   return retval;
  325. }
  326.  
  327. Matrix
  328. Matrix::stack (const Matrix& a) const
  329. {
  330.   int nr = rows ();
  331.   int nc = cols ();
  332.   if (nc != a.cols ())
  333.     {
  334.       (*current_liboctave_error_handler)
  335.     ("column dimension mismatch for stack");
  336.       return Matrix ();
  337.     }
  338.  
  339.   int nr_insert = nr;
  340.   Matrix retval (nr + a.rows (), nc);
  341.   retval.insert (*this, 0, 0);
  342.   retval.insert (a, nr_insert, 0);
  343.   return retval;
  344. }
  345.  
  346. Matrix
  347. Matrix::stack (const RowVector& a) const
  348. {
  349.   int nr = rows ();
  350.   int nc = cols ();
  351.   if (nc != a.length ())
  352.     {
  353.       (*current_liboctave_error_handler)
  354.     ("column dimension mismatch for stack");
  355.       return Matrix ();
  356.     }
  357.  
  358.   int nr_insert = nr;
  359.   Matrix retval (nr + 1, nc);
  360.   retval.insert (*this, 0, 0);
  361.   retval.insert (a, nr_insert, 0);
  362.   return retval;
  363. }
  364.  
  365. Matrix
  366. Matrix::stack (const ColumnVector& a) const
  367. {
  368.   int nr = rows ();
  369.   int nc = cols ();
  370.   if (nc != 1)
  371.     {
  372.       (*current_liboctave_error_handler)
  373.     ("column dimension mismatch for stack");
  374.       return Matrix ();
  375.     }
  376.  
  377.   int nr_insert = nr;
  378.   Matrix retval (nr + a.length (), nc);
  379.   retval.insert (*this, 0, 0);
  380.   retval.insert (a, nr_insert, 0);
  381.   return retval;
  382. }
  383.  
  384. Matrix
  385. Matrix::stack (const DiagMatrix& a) const
  386. {
  387.   int nr = rows ();
  388.   int nc = cols ();
  389.   if (nc != a.cols ())
  390.     {
  391.       (*current_liboctave_error_handler)
  392.     ("column dimension mismatch for stack");
  393.       return Matrix ();
  394.     }
  395.  
  396.   int nr_insert = nr;
  397.   Matrix retval (nr + a.rows (), nc);
  398.   retval.insert (*this, 0, 0);
  399.   retval.insert (a, nr_insert, 0);
  400.   return retval;
  401. }
  402.  
  403. Matrix
  404. Matrix::transpose (void) const
  405. {
  406.   int nr = rows ();
  407.   int nc = cols ();
  408.   Matrix result (nc, nr);
  409.   if (length () > 0)
  410.     {
  411.       for (int j = 0; j < nc; j++)
  412.     for (int i = 0; i < nr; i++)
  413.       result.elem (j, i) = elem (i, j);
  414.     }
  415.   return result;
  416. }
  417.  
  418. Matrix
  419. real (const ComplexMatrix& a)
  420. {
  421.   int a_len = a.length ();
  422.   Matrix retval;
  423.   if (a_len > 0)
  424.     retval = Matrix (real_dup (a.data (), a_len), a.rows (), a.cols ());
  425.   return retval;
  426. }
  427.  
  428. Matrix
  429. imag (const ComplexMatrix& a)
  430. {
  431.   int a_len = a.length ();
  432.   Matrix retval;
  433.   if (a_len > 0)
  434.     retval = Matrix (imag_dup (a.data (), a_len), a.rows (), a.cols ());
  435.   return retval;
  436. }
  437.  
  438. Matrix
  439. Matrix::extract (int r1, int c1, int r2, int c2) const
  440. {
  441.   if (r1 > r2) { int tmp = r1; r1 = r2; r2 = tmp; }
  442.   if (c1 > c2) { int tmp = c1; c1 = c2; c2 = tmp; }
  443.  
  444.   int new_r = r2 - r1 + 1;
  445.   int new_c = c2 - c1 + 1;
  446.  
  447.   Matrix result (new_r, new_c);
  448.  
  449.   for (int j = 0; j < new_c; j++)
  450.     for (int i = 0; i < new_r; i++)
  451.       result.elem (i, j) = elem (r1+i, c1+j);
  452.  
  453.   return result;
  454. }
  455.  
  456. // extract row or column i.
  457.  
  458. RowVector
  459. Matrix::row (int i) const
  460. {
  461.   int nc = cols ();
  462.   if (i < 0 || i >= rows ())
  463.     {
  464.       (*current_liboctave_error_handler) ("invalid row selection");
  465.       return RowVector ();
  466.     }
  467.  
  468.   RowVector retval (nc);
  469.   for (int j = 0; j < nc; j++)
  470.     retval.elem (j) = elem (i, j);
  471.  
  472.   return retval;
  473. }
  474.  
  475. RowVector
  476. Matrix::row (char *s) const
  477. {
  478.   if (! s)
  479.     {
  480.       (*current_liboctave_error_handler) ("invalid row selection");
  481.       return RowVector ();
  482.     }
  483.  
  484.   char c = *s;
  485.   if (c == 'f' || c == 'F')
  486.     return row (0);
  487.   else if (c == 'l' || c == 'L')
  488.     return row (rows () - 1);
  489.   else
  490.     {
  491.       (*current_liboctave_error_handler) ("invalid row selection");
  492.       return RowVector ();
  493.     }
  494. }
  495.  
  496. ColumnVector
  497. Matrix::column (int i) const
  498. {
  499.   int nr = rows ();
  500.   if (i < 0 || i >= cols ())
  501.     {
  502.       (*current_liboctave_error_handler) ("invalid column selection");
  503.       return ColumnVector ();
  504.     }
  505.  
  506.   ColumnVector retval (nr);
  507.   for (int j = 0; j < nr; j++)
  508.     retval.elem (j) = elem (j, i);
  509.  
  510.   return retval;
  511. }
  512.  
  513. ColumnVector
  514. Matrix::column (char *s) const
  515. {
  516.   if (! s)
  517.     {
  518.       (*current_liboctave_error_handler) ("invalid column selection");
  519.       return ColumnVector ();
  520.     }
  521.  
  522.   char c = *s;
  523.   if (c == 'f' || c == 'F')
  524.     return column (0);
  525.   else if (c == 'l' || c == 'L')
  526.     return column (cols () - 1);
  527.   else
  528.     {
  529.       (*current_liboctave_error_handler) ("invalid column selection");
  530.       return ColumnVector ();
  531.     }
  532. }
  533.  
  534. Matrix
  535. Matrix::inverse (void) const
  536. {
  537.   int info;
  538.   double rcond;
  539.   return inverse (info, rcond);
  540. }
  541.  
  542. Matrix
  543. Matrix::inverse (int& info) const
  544. {
  545.   double rcond;
  546.   return inverse (info, rcond);
  547. }
  548.  
  549. Matrix
  550. Matrix::inverse (int& info, double& rcond, int force) const
  551. {
  552.   Matrix retval;
  553.  
  554.   int nr = rows ();
  555.   int nc = cols ();
  556.  
  557.   if (nr != nc || nr == 0 || nc == 0)
  558.     (*current_liboctave_error_handler) ("inverse requires square matrix");
  559.   else
  560.     {
  561.       info = 0;
  562.  
  563.       Array<int> ipvt (nr);
  564.       int *pipvt = ipvt.fortran_vec ();
  565.  
  566.       Array<double> z (nr);
  567.       double *pz = z.fortran_vec ();
  568.  
  569.       retval = *this;
  570.       double *tmp_data = retval.fortran_vec ();
  571.  
  572.       F77_XFCN (dgeco, DGECO, (tmp_data, nr, nc, pipvt, rcond, pz));
  573.  
  574.       if (f77_exception_encountered)
  575.     (*current_liboctave_error_handler) ("unrecoverable error in dgeco");
  576.       else
  577.     {
  578.       volatile double rcond_plus_one = rcond + 1.0;
  579.  
  580.       if (rcond_plus_one == 1.0)
  581.         info = -1;
  582.  
  583.       if (info == -1 && ! force)
  584.         retval = *this; // Restore matrix contents.
  585.       else
  586.         {
  587.           double *dummy = 0;
  588.  
  589.           F77_XFCN (dgedi, DGEDI, (tmp_data, nr, nc, pipvt, dummy,
  590.                        pz, 1));
  591.  
  592.           if (f77_exception_encountered)
  593.         (*current_liboctave_error_handler)
  594.           ("unrecoverable error in dgedi");
  595.         }
  596.     }
  597.     }
  598.  
  599.   return retval;
  600. }
  601.  
  602. Matrix
  603. Matrix::pseudo_inverse (double tol)
  604. {
  605.   SVD result (*this);
  606.  
  607.   DiagMatrix S = result.singular_values ();
  608.   Matrix U = result.left_singular_matrix ();
  609.   Matrix V = result.right_singular_matrix ();
  610.  
  611.   ColumnVector sigma = S.diag ();
  612.  
  613.   int r = sigma.length () - 1;
  614.   int nr = rows ();
  615.   int nc = cols ();
  616.  
  617.   if (tol <= 0.0)
  618.     {
  619.       if (nr > nc)
  620.     tol = nr * sigma.elem (0) * DBL_EPSILON;
  621.       else
  622.     tol = nc * sigma.elem (0) * DBL_EPSILON;
  623.     }
  624.  
  625.   while (r >= 0 && sigma.elem (r) < tol)
  626.     r--;
  627.  
  628.   if (r < 0)
  629.     return Matrix (nc, nr, 0.0);
  630.   else
  631.     {
  632.       Matrix Ur = U.extract (0, 0, nr-1, r);
  633.       DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
  634.       Matrix Vr = V.extract (0, 0, nc-1, r);
  635.       return Vr * D * Ur.transpose ();
  636.     }
  637. }
  638.  
  639. ComplexMatrix
  640. Matrix::fourier (void) const
  641. {
  642.   ComplexMatrix retval;
  643.  
  644.   int nr = rows ();
  645.   int nc = cols ();
  646.  
  647.   int npts, nsamples;
  648.  
  649.   if (nr == 1 || nc == 1)
  650.     {
  651.       npts = nr > nc ? nr : nc;
  652.       nsamples = 1;
  653.     }
  654.   else
  655.     {
  656.       npts = nr;
  657.       nsamples = nc;
  658.     }
  659.  
  660.   int nn = 4*npts+15;
  661.  
  662.   Array<Complex> wsave (nn);
  663.   Complex *pwsave = wsave.fortran_vec ();
  664.  
  665.   retval = *this;
  666.   Complex *tmp_data = retval.fortran_vec ();
  667.  
  668.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  669.  
  670.   for (int j = 0; j < nsamples; j++)
  671.     F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
  672.  
  673.   return retval;
  674. }
  675.  
  676. ComplexMatrix
  677. Matrix::ifourier (void) const
  678. {
  679.   ComplexMatrix retval;
  680.  
  681.   int nr = rows ();
  682.   int nc = cols ();
  683.  
  684.   int npts, nsamples;
  685.  
  686.   if (nr == 1 || nc == 1)
  687.     {
  688.       npts = nr > nc ? nr : nc;
  689.       nsamples = 1;
  690.     }
  691.   else
  692.     {
  693.       npts = nr;
  694.       nsamples = nc;
  695.     }
  696.  
  697.   int nn = 4*npts+15;
  698.  
  699.   Array<Complex> wsave (nn);
  700.   Complex *pwsave = wsave.fortran_vec ();
  701.  
  702.   retval = *this;
  703.   Complex *tmp_data = retval.fortran_vec ();
  704.  
  705.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  706.  
  707.   for (int j = 0; j < nsamples; j++)
  708.     F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
  709.  
  710.   for (int j = 0; j < npts*nsamples; j++)
  711.     tmp_data[j] = tmp_data[j] / (double) npts;
  712.  
  713.   return retval;
  714. }
  715.  
  716. ComplexMatrix
  717. Matrix::fourier2d (void) const
  718. {
  719.   ComplexMatrix retval;
  720.  
  721.   int nr = rows ();
  722.   int nc = cols ();
  723.  
  724.   int npts, nsamples;
  725.  
  726.   if (nr == 1 || nc == 1)
  727.     {
  728.       npts = nr > nc ? nr : nc;
  729.       nsamples = 1;
  730.     }
  731.   else
  732.     {
  733.       npts = nr;
  734.       nsamples = nc;
  735.     }
  736.  
  737.   int nn = 4*npts+15;
  738.  
  739.   Array<Complex> wsave (nn);
  740.   Complex *pwsave = wsave.fortran_vec ();
  741.  
  742.   retval = *this;
  743.   Complex *tmp_data = retval.fortran_vec ();
  744.  
  745.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  746.  
  747.   for (int j = 0; j < nsamples; j++)
  748.     F77_FCN (cfftf, CFFTF) (npts, &tmp_data[npts*j], pwsave);
  749.  
  750.   npts = nc;
  751.   nsamples = nr;
  752.   nn = 4*npts+15;
  753.  
  754.   wsave.resize (nn);
  755.   pwsave = wsave.fortran_vec ();
  756.  
  757.   Array<Complex> row (npts);
  758.   Complex *prow = row.fortran_vec ();
  759.  
  760.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  761.  
  762.   for (int j = 0; j < nsamples; j++)
  763.     {
  764.       for (int i = 0; i < npts; i++)
  765.     prow[i] = tmp_data[i*nr + j];
  766.  
  767.       F77_FCN (cfftf, CFFTF) (npts, prow, pwsave);
  768.  
  769.       for (int i = 0; i < npts; i++)
  770.     tmp_data[i*nr + j] = prow[i];
  771.     }
  772.  
  773.   return retval;
  774. }
  775.  
  776. ComplexMatrix
  777. Matrix::ifourier2d (void) const
  778. {
  779.   ComplexMatrix retval;
  780.  
  781.   int nr = rows ();
  782.   int nc = cols ();
  783.  
  784.   int npts, nsamples;
  785.  
  786.   if (nr == 1 || nc == 1)
  787.     {
  788.       npts = nr > nc ? nr : nc;
  789.       nsamples = 1;
  790.     }
  791.   else
  792.     {
  793.       npts = nr;
  794.       nsamples = nc;
  795.     }
  796.  
  797.   int nn = 4*npts+15;
  798.  
  799.   Array<Complex> wsave (nn);
  800.   Complex *pwsave = wsave.fortran_vec ();
  801.  
  802.   retval = *this;
  803.   Complex *tmp_data = retval.fortran_vec ();
  804.  
  805.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  806.  
  807.   for (int j = 0; j < nsamples; j++)
  808.     F77_FCN (cfftb, CFFTB) (npts, &tmp_data[npts*j], pwsave);
  809.  
  810.   for (int j = 0; j < npts*nsamples; j++)
  811.     tmp_data[j] = tmp_data[j] / (double) npts;
  812.  
  813.   npts = nc;
  814.   nsamples = nr;
  815.   nn = 4*npts+15;
  816.  
  817.   wsave.resize (nn);
  818.   pwsave = wsave.fortran_vec ();
  819.  
  820.   Array<Complex> row (npts);
  821.   Complex *prow = row.fortran_vec ();
  822.  
  823.   F77_FCN (cffti, CFFTI) (npts, pwsave);
  824.  
  825.   for (int j = 0; j < nsamples; j++)
  826.     {
  827.       for (int i = 0; i < npts; i++)
  828.     prow[i] = tmp_data[i*nr + j];
  829.  
  830.       F77_FCN (cfftb, CFFTB) (npts, prow, pwsave);
  831.  
  832.       for (int i = 0; i < npts; i++)
  833.     tmp_data[i*nr + j] = prow[i] / (double) npts;
  834.     }
  835.  
  836.   return retval;
  837. }
  838.  
  839. DET
  840. Matrix::determinant (void) const
  841. {
  842.   int info;
  843.   double rcond;
  844.   return determinant (info, rcond);
  845. }
  846.  
  847. DET
  848. Matrix::determinant (int& info) const
  849. {
  850.   double rcond;
  851.   return determinant (info, rcond);
  852. }
  853.  
  854. DET
  855. Matrix::determinant (int& info, double& rcond) const
  856. {
  857.   DET retval;
  858.  
  859.   int nr = rows ();
  860.   int nc = cols ();
  861.  
  862.   if (nr == 0 || nc == 0)
  863.     {
  864.       double d[2];
  865.       d[0] = 1.0;
  866.       d[1] = 0.0;
  867.       retval = DET (d);
  868.     }
  869.   else
  870.     {
  871.       info = 0;
  872.  
  873.       Array<int> ipvt (nr);
  874.       int *pipvt = ipvt.fortran_vec ();
  875.  
  876.       Array<double> z (nr);
  877.       double *pz = z.fortran_vec ();
  878.  
  879.       Matrix atmp = *this;
  880.       double *tmp_data = atmp.fortran_vec ();
  881.  
  882.       F77_XFCN (dgeco, DGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
  883.  
  884.       if (f77_exception_encountered)
  885.     (*current_liboctave_error_handler) ("unrecoverable error in dgeco");
  886.       else
  887.     {
  888.       volatile double rcond_plus_one = rcond + 1.0;
  889.  
  890.       if (rcond_plus_one == 1.0)
  891.         {
  892.           info = -1;
  893.           retval = DET ();
  894.         }
  895.       else
  896.         {
  897.           double d[2];
  898.  
  899.           F77_XFCN (dgedi, DGEDI, (tmp_data, nr, nr, pipvt, d, pz, 10));
  900.  
  901.           if (f77_exception_encountered)
  902.         (*current_liboctave_error_handler)
  903.           ("unrecoverable error in dgedi");
  904.           else
  905.         retval = DET (d);
  906.         }
  907.     }
  908.     }
  909.  
  910.   return retval;
  911. }
  912.  
  913. Matrix
  914. Matrix::solve (const Matrix& b) const
  915. {
  916.   int info;
  917.   double rcond;
  918.   return solve (b, info, rcond);
  919. }
  920.  
  921. Matrix
  922. Matrix::solve (const Matrix& b, int& info) const
  923. {
  924.   double rcond;
  925.   return solve (b, info, rcond);
  926. }
  927.  
  928. Matrix
  929. Matrix::solve (const Matrix& b, int& info, double& rcond) const
  930. {
  931.   Matrix retval;
  932.  
  933.   int nr = rows ();
  934.   int nc = cols ();
  935.  
  936.   if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
  937.     (*current_liboctave_error_handler)
  938.       ("matrix dimension mismatch solution of linear equations");
  939.   else
  940.     {
  941.       info = 0;
  942.  
  943.       Array<int> ipvt (nr);
  944.       int *pipvt = ipvt.fortran_vec ();
  945.  
  946.       Array<double> z (nr);
  947.       double *pz = z.fortran_vec ();
  948.  
  949.       Matrix atmp = *this;
  950.       double *tmp_data = atmp.fortran_vec ();
  951.  
  952.       F77_XFCN (dgeco, DGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
  953.  
  954.       if (f77_exception_encountered)
  955.     (*current_liboctave_error_handler) ("unrecoverable error in dgeco");
  956.       else
  957.     {
  958.       volatile double rcond_plus_one = rcond + 1.0;
  959.  
  960.       if (rcond_plus_one == 1.0)
  961.         {
  962.           info = -2;
  963.         }
  964.       else
  965.         {
  966.           retval = b;
  967.           double *result = retval.fortran_vec ();
  968.  
  969.           int b_nc = b.cols ();
  970.  
  971.           for (volatile int j = 0; j < b_nc; j++)
  972.         {
  973.           F77_XFCN (dgesl, DGESL, (tmp_data, nr, nr, pipvt,
  974.                        &result[nr*j], 0)); 
  975.  
  976.           if (f77_exception_encountered)
  977.             {
  978.               (*current_liboctave_error_handler)
  979.             ("unrecoverable error in dgesl");
  980.  
  981.               break;
  982.             }
  983.         }
  984.         }
  985.     }
  986.     }
  987.  
  988.   return retval;
  989. }
  990.  
  991. ComplexMatrix
  992. Matrix::solve (const ComplexMatrix& b) const
  993. {
  994.   ComplexMatrix tmp (*this);
  995.   return tmp.solve (b);
  996. }
  997.  
  998. ComplexMatrix
  999. Matrix::solve (const ComplexMatrix& b, int& info) const
  1000. {
  1001.   ComplexMatrix tmp (*this);
  1002.   return tmp.solve (b, info);
  1003. }
  1004.  
  1005. ComplexMatrix
  1006. Matrix::solve (const ComplexMatrix& b, int& info, double& rcond) const
  1007. {
  1008.   ComplexMatrix tmp (*this);
  1009.   return tmp.solve (b, info, rcond);
  1010. }
  1011.  
  1012. ColumnVector
  1013. Matrix::solve (const ColumnVector& b) const
  1014. {
  1015.   int info; double rcond;
  1016.   return solve (b, info, rcond);
  1017. }
  1018.  
  1019. ColumnVector
  1020. Matrix::solve (const ColumnVector& b, int& info) const
  1021. {
  1022.   double rcond;
  1023.   return solve (b, info, rcond);
  1024. }
  1025.  
  1026. ColumnVector
  1027. Matrix::solve (const ColumnVector& b, int& info, double& rcond) const
  1028. {
  1029.   ColumnVector retval;
  1030.  
  1031.   int nr = rows ();
  1032.   int nc = cols ();
  1033.  
  1034.   if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
  1035.     (*current_liboctave_error_handler)
  1036.       ("matrix dimension mismatch solution of linear equations");
  1037.   else
  1038.     {
  1039.       info = 0;
  1040.  
  1041.       Array<int> ipvt (nr);
  1042.       int *pipvt = ipvt.fortran_vec ();
  1043.  
  1044.       Array<double> z (nr);
  1045.       double *pz = z.fortran_vec ();
  1046.  
  1047.       Matrix atmp = *this;
  1048.       double *tmp_data = atmp.fortran_vec ();
  1049.  
  1050.       F77_XFCN (dgeco, DGECO, (tmp_data, nr, nr, pipvt, rcond, pz));
  1051.  
  1052.       if (f77_exception_encountered)
  1053.     (*current_liboctave_error_handler)
  1054.       ("unrecoverable error in dgeco");
  1055.       else
  1056.     {
  1057.       volatile double rcond_plus_one = rcond + 1.0;
  1058.  
  1059.       if (rcond_plus_one == 1.0)
  1060.         {
  1061.           info = -2;
  1062.         }
  1063.       else
  1064.         {
  1065.           retval = b;
  1066.           double *result = retval.fortran_vec ();
  1067.  
  1068.           F77_XFCN (dgesl, DGESL, (tmp_data, nr, nr, pipvt, result, 0));
  1069.  
  1070.           if (f77_exception_encountered)
  1071.         (*current_liboctave_error_handler)
  1072.           ("unrecoverable error in dgesl");
  1073.         }
  1074.     }
  1075.     }
  1076.  
  1077.   return retval;
  1078. }
  1079.  
  1080. ComplexColumnVector
  1081. Matrix::solve (const ComplexColumnVector& b) const
  1082. {
  1083.   ComplexMatrix tmp (*this);
  1084.   return tmp.solve (b);
  1085. }
  1086.  
  1087. ComplexColumnVector
  1088. Matrix::solve (const ComplexColumnVector& b, int& info) const
  1089. {
  1090.   ComplexMatrix tmp (*this);
  1091.   return tmp.solve (b, info);
  1092. }
  1093.  
  1094. ComplexColumnVector
  1095. Matrix::solve (const ComplexColumnVector& b, int& info, double& rcond) const
  1096. {
  1097.   ComplexMatrix tmp (*this);
  1098.   return tmp.solve (b, info, rcond);
  1099. }
  1100.  
  1101. Matrix
  1102. Matrix::lssolve (const Matrix& b) const
  1103. {
  1104.   int info;
  1105.   int rank;
  1106.   return lssolve (b, info, rank);
  1107. }
  1108.  
  1109. Matrix
  1110. Matrix::lssolve (const Matrix& b, int& info) const
  1111. {
  1112.   int rank;
  1113.   return lssolve (b, info, rank);
  1114. }
  1115.  
  1116. Matrix
  1117. Matrix::lssolve (const Matrix& b, int& info, int& rank) const
  1118. {
  1119.   Matrix retval;
  1120.  
  1121.   int nrhs = b.cols ();
  1122.  
  1123.   int m = rows ();
  1124.   int n = cols ();
  1125.  
  1126.   if (m == 0 || n == 0 || m != b.rows ())
  1127.     (*current_liboctave_error_handler)
  1128.       ("matrix dimension mismatch in solution of least squares problem");
  1129.   else
  1130.     {
  1131.       Matrix atmp = *this;
  1132.       double *tmp_data = atmp.fortran_vec ();
  1133.  
  1134.       int nrr = m > n ? m : n;
  1135.       Matrix result (nrr, nrhs);
  1136.  
  1137.       for (int j = 0; j < nrhs; j++)
  1138.     for (int i = 0; i < m; i++)
  1139.       result.elem (i, j) = b.elem (i, j);
  1140.  
  1141.       double *presult = result.fortran_vec ();
  1142.  
  1143.       int len_s = m < n ? m : n;
  1144.       Array<double> s (len_s);
  1145.       double *ps = s.fortran_vec ();
  1146.  
  1147.       double rcond = -1.0;
  1148.  
  1149.       int lwork;
  1150.       if (m < n)
  1151.     lwork = 3*m + (2*m > nrhs
  1152.                ? (2*m > n ? 2*m : n)
  1153.                : (nrhs > n ? nrhs : n));
  1154.       else
  1155.     lwork = 3*n + (2*n > nrhs
  1156.                ? (2*n > m ? 2*n : m)
  1157.                : (nrhs > m ? nrhs : m));
  1158.  
  1159.       lwork *= 16;
  1160.  
  1161.       Array<double> work (lwork);
  1162.       double *pwork = work.fortran_vec ();
  1163.  
  1164.       F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult, nrr, ps,
  1165.                  rcond, rank, pwork, lwork, info));
  1166.  
  1167.       if (f77_exception_encountered)
  1168.     (*current_liboctave_error_handler) ("unrecoverable error in dgelss");
  1169.       else
  1170.     {
  1171.       retval.resize (n, nrhs);
  1172.       for (int j = 0; j < nrhs; j++)
  1173.         for (int i = 0; i < n; i++)
  1174.           retval.elem (i, j) = result.elem (i, j);
  1175.     }
  1176.     }
  1177.  
  1178.   return retval;
  1179. }
  1180.  
  1181. ComplexMatrix
  1182. Matrix::lssolve (const ComplexMatrix& b) const
  1183. {
  1184.   ComplexMatrix tmp (*this);
  1185.   int info;
  1186.   int rank;
  1187.   return tmp.lssolve (b, info, rank);
  1188. }
  1189.  
  1190. ComplexMatrix
  1191. Matrix::lssolve (const ComplexMatrix& b, int& info) const
  1192. {
  1193.   ComplexMatrix tmp (*this);
  1194.   int rank;
  1195.   return tmp.lssolve (b, info, rank);
  1196. }
  1197.  
  1198. ComplexMatrix
  1199. Matrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const
  1200. {
  1201.   ComplexMatrix tmp (*this);
  1202.   return tmp.lssolve (b, info, rank);
  1203. }
  1204.  
  1205. ColumnVector
  1206. Matrix::lssolve (const ColumnVector& b) const
  1207. {
  1208.   int info;
  1209.   int rank;
  1210.   return lssolve (b, info, rank);
  1211. }
  1212.  
  1213. ColumnVector
  1214. Matrix::lssolve (const ColumnVector& b, int& info) const
  1215. {
  1216.   int rank;
  1217.   return lssolve (b, info, rank);
  1218. }
  1219.  
  1220. ColumnVector
  1221. Matrix::lssolve (const ColumnVector& b, int& info, int& rank) const
  1222. {
  1223.   ColumnVector retval;
  1224.  
  1225.   int nrhs = 1;
  1226.  
  1227.   int m = rows ();
  1228.   int n = cols ();
  1229.  
  1230.   if (m == 0 || n == 0 || m != b.length ())
  1231.     (*current_liboctave_error_handler)
  1232.       ("matrix dimension mismatch in solution of least squares problem");
  1233.   else
  1234.     {
  1235.       Matrix atmp = *this;
  1236.       double *tmp_data = atmp.fortran_vec ();
  1237.  
  1238.       int nrr = m > n ? m : n;
  1239.       ColumnVector result (nrr);
  1240.  
  1241.       for (int i = 0; i < m; i++)
  1242.     result.elem (i) = b.elem (i);
  1243.  
  1244.       double *presult = result.fortran_vec ();
  1245.  
  1246.       int len_s = m < n ? m : n;
  1247.       Array<double> s (len_s);
  1248.       double *ps = s.fortran_vec ();
  1249.  
  1250.       double rcond = -1.0;
  1251.  
  1252.       int lwork;
  1253.       if (m < n)
  1254.     lwork = 3*m + (2*m > nrhs
  1255.                ? (2*m > n ? 2*m : n)
  1256.                : (nrhs > n ? nrhs : n));
  1257.       else
  1258.     lwork = 3*n + (2*n > nrhs
  1259.                ? (2*n > m ? 2*n : m)
  1260.                : (nrhs > m ? nrhs : m));
  1261.  
  1262.       lwork *= 16;
  1263.  
  1264.       Array<double> work (lwork);
  1265.       double *pwork = work.fortran_vec ();
  1266.  
  1267.       F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult, nrr,
  1268.                  ps, rcond, rank, pwork, lwork, info));
  1269.  
  1270.       if (f77_exception_encountered)
  1271.     (*current_liboctave_error_handler) ("unrecoverable error in dgelss");
  1272.       else
  1273.     {
  1274.       retval.resize (n);
  1275.       for (int i = 0; i < n; i++)
  1276.         retval.elem (i) = result.elem (i);
  1277.     }
  1278.     }
  1279.  
  1280.   return retval;
  1281. }
  1282.  
  1283. ComplexColumnVector
  1284. Matrix::lssolve (const ComplexColumnVector& b) const
  1285. {
  1286.   ComplexMatrix tmp (*this);
  1287.   return tmp.lssolve (b);
  1288. }
  1289.  
  1290. ComplexColumnVector
  1291. Matrix::lssolve (const ComplexColumnVector& b, int& info) const
  1292. {
  1293.   ComplexMatrix tmp (*this);
  1294.   return tmp.lssolve (b, info);
  1295. }
  1296.  
  1297. ComplexColumnVector
  1298. Matrix::lssolve (const ComplexColumnVector& b, int& info, int& rank) const
  1299. {
  1300.   ComplexMatrix tmp (*this);
  1301.   return tmp.lssolve (b, info, rank);
  1302. }
  1303.  
  1304. // Constants for matrix exponential calculation.
  1305.  
  1306. static double padec [] =
  1307. {
  1308.   5.0000000000000000e-1,
  1309.   1.1666666666666667e-1,
  1310.   1.6666666666666667e-2,
  1311.   1.6025641025641026e-3,
  1312.   1.0683760683760684e-4,
  1313.   4.8562548562548563e-6,
  1314.   1.3875013875013875e-7,
  1315.   1.9270852604185938e-9,
  1316. };
  1317.  
  1318. Matrix
  1319. Matrix::expm (void) const
  1320. {
  1321.   Matrix retval;
  1322.  
  1323.   Matrix m = *this;
  1324.  
  1325.   int nc = columns ();
  1326.  
  1327.   // trace shift value
  1328.   double trshift = 0;
  1329.  
  1330.   // Preconditioning step 1: trace normalization.
  1331.  
  1332.   for (int i = 0; i < nc; i++)
  1333.     trshift += m.elem (i, i);
  1334.  
  1335.   trshift /= nc;
  1336.  
  1337.   for (int i = 0; i < nc; i++)
  1338.     m.elem (i, i) -= trshift;
  1339.  
  1340.   // Preconditioning step 2: balancing.
  1341.  
  1342.   AEPBALANCE mbal (m, "B");
  1343.   m = mbal.balanced_matrix ();
  1344.   Matrix d = mbal.balancing_matrix ();
  1345.  
  1346.   // Preconditioning step 3: scaling.
  1347.  
  1348.   ColumnVector work(nc);
  1349.   double inf_norm
  1350.     = F77_FCN (dlange, DLANGE) ("I", nc, nc, m.fortran_vec (),nc,
  1351.                 work.fortran_vec ());
  1352.  
  1353.   int sqpow = (int) (inf_norm > 0.0
  1354.              ? (1.0 + log (inf_norm) / log (2.0))
  1355.              : 0.0);
  1356.  
  1357.   // Check whether we need to square at all.
  1358.  
  1359.   if (sqpow < 0)
  1360.     sqpow = 0;
  1361.  
  1362.   if (sqpow > 0)
  1363.     {
  1364.       double scale_factor = 1.0;
  1365.       for (int i = 0; i < sqpow; i++)
  1366.     scale_factor *= 2.0;
  1367.  
  1368.       m = m / scale_factor;
  1369.     }
  1370.  
  1371.   // npp, dpp: pade' approx polynomial matrices.
  1372.  
  1373.   Matrix npp (nc, nc, 0.0);
  1374.   Matrix dpp = npp;
  1375.  
  1376.   // Now powers a^8 ... a^1.
  1377.  
  1378.   int minus_one_j = -1;
  1379.   for (int j = 7; j >= 0; j--)
  1380.     {
  1381.       npp = m * npp + m * padec[j];
  1382.       dpp = m * dpp + m * (minus_one_j * padec[j]);
  1383.       minus_one_j *= -1;
  1384.     }
  1385.  
  1386.   // Zero power.
  1387.  
  1388.   dpp = -dpp;
  1389.   for(int j = 0; j < nc; j++)
  1390.     {
  1391.       npp.elem (j, j) += 1.0;
  1392.       dpp.elem (j, j) += 1.0;
  1393.     }
  1394.  
  1395.   // Compute pade approximation = inverse (dpp) * npp.
  1396.  
  1397.   retval = dpp.solve (npp);
  1398.  
  1399.   // Reverse preconditioning step 3: repeated squaring.
  1400.  
  1401.   while (sqpow)
  1402.     {
  1403.       retval = retval * retval;
  1404.       sqpow--;
  1405.     }
  1406.  
  1407.   // Reverse preconditioning step 2: inverse balancing.
  1408.  
  1409.   retval = retval.transpose();
  1410.   d = d.transpose ();
  1411.   retval = retval * d;
  1412.   retval = d.solve (retval);
  1413.   retval = retval.transpose ();
  1414.  
  1415.   // Reverse preconditioning step 1: fix trace normalization.
  1416.  
  1417.   return retval * exp (trshift);
  1418. }
  1419.  
  1420. Matrix&
  1421. Matrix::operator += (const Matrix& a)
  1422. {
  1423.   int nr = rows ();
  1424.   int nc = cols ();
  1425.  
  1426.   int a_nr = a.rows ();
  1427.   int a_nc = a.cols ();
  1428.  
  1429.   if (nr != a_nr || nc != a_nc)
  1430.     {
  1431.       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
  1432.       return *this;
  1433.     }
  1434.  
  1435.   if (nr == 0 || nc == 0)
  1436.     return *this;
  1437.  
  1438.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  1439.  
  1440.   add2 (d, a.data (), length ());
  1441.  
  1442.   return *this;
  1443. }
  1444.  
  1445. Matrix&
  1446. Matrix::operator -= (const Matrix& a)
  1447. {
  1448.   int nr = rows ();
  1449.   int nc = cols ();
  1450.  
  1451.   int a_nr = a.rows ();
  1452.   int a_nc = a.cols ();
  1453.  
  1454.   if (nr != a_nr || nc != a_nc)
  1455.     {
  1456.       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
  1457.       return *this;
  1458.     }
  1459.  
  1460.   if (nr == 0 || nc == 0)
  1461.     return *this;
  1462.  
  1463.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  1464.  
  1465.   subtract2 (d, a.data (), length ());
  1466.  
  1467.   return *this;
  1468. }
  1469.  
  1470. Matrix&
  1471. Matrix::operator += (const DiagMatrix& a)
  1472. {
  1473.   int nr = rows ();
  1474.   int nc = cols ();
  1475.  
  1476.   int a_nr = a.rows ();
  1477.   int a_nc = a.cols ();
  1478.  
  1479.   if (nr != a_nr || nc != a_nc)
  1480.     {
  1481.       gripe_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
  1482.       return *this;
  1483.     }
  1484.  
  1485.   for (int i = 0; i < a.length (); i++)
  1486.     elem (i, i) += a.elem (i, i);
  1487.  
  1488.   return *this;
  1489. }
  1490.  
  1491. Matrix&
  1492. Matrix::operator -= (const DiagMatrix& a)
  1493. {
  1494.   int nr = rows ();
  1495.   int nc = cols ();
  1496.  
  1497.   int a_nr = a.rows ();
  1498.   int a_nc = a.cols ();
  1499.  
  1500.   if (nr != a_nr || nc != a_nc)
  1501.     {
  1502.       gripe_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
  1503.       return *this;
  1504.     }
  1505.  
  1506.   for (int i = 0; i < a.length (); i++)
  1507.     elem (i, i) -= a.elem (i, i);
  1508.  
  1509.   return *this;
  1510. }
  1511.  
  1512. // unary operations
  1513.  
  1514. Matrix
  1515. Matrix::operator ! (void) const
  1516. {
  1517.   int nr = rows ();
  1518.   int nc = cols ();
  1519.  
  1520.   Matrix b (nr, nc);
  1521.  
  1522.   for (int j = 0; j < nc; j++)
  1523.     for (int i = 0; i < nr; i++)
  1524.       b.elem (i, j) = ! elem (i, j);
  1525.  
  1526.   return b;
  1527. }
  1528.  
  1529. // column vector by row vector -> matrix operations
  1530.  
  1531. Matrix
  1532. operator * (const ColumnVector& v, const RowVector& a)
  1533. {
  1534.   Matrix retval;
  1535.  
  1536.   int len = v.length ();
  1537.   int a_len = a.length ();
  1538.  
  1539.   if (len != a_len)
  1540.     gripe_nonconformant ("operator *", len, 1, 1, a_len);
  1541.   else
  1542.     {
  1543.       if (len != 0)
  1544.     {
  1545.       retval.resize (len, a_len);
  1546.       double *c = retval.fortran_vec ();
  1547.  
  1548.       F77_XFCN (dgemm, DGEMM, ("N", "N", len, a_len, 1, 1.0,
  1549.                    v.data (), len, a.data (), 1, 0.0,
  1550.                    c, len, 1L, 1L));
  1551.  
  1552.       if (f77_exception_encountered)
  1553.         (*current_liboctave_error_handler)
  1554.           ("unrecoverable error in dgemm");
  1555.     }
  1556.     }
  1557.  
  1558.   return retval;
  1559. }
  1560.  
  1561. // diagonal matrix by scalar -> matrix operations
  1562.  
  1563. Matrix
  1564. operator + (const DiagMatrix& a, double s)
  1565. {
  1566.   Matrix tmp (a.rows (), a.cols (), s);
  1567.   return a + tmp;
  1568. }
  1569.  
  1570. Matrix
  1571. operator - (const DiagMatrix& a, double s)
  1572. {
  1573.   Matrix tmp (a.rows (), a.cols (), -s);
  1574.   return a + tmp;
  1575. }
  1576.  
  1577. // scalar by diagonal matrix -> matrix operations
  1578.  
  1579. Matrix
  1580. operator + (double s, const DiagMatrix& a)
  1581. {
  1582.   Matrix tmp (a.rows (), a.cols (), s);
  1583.   return tmp + a;
  1584. }
  1585.  
  1586. Matrix
  1587. operator - (double s, const DiagMatrix& a)
  1588. {
  1589.   Matrix tmp (a.rows (), a.cols (), s);
  1590.   return tmp - a;
  1591. }
  1592.  
  1593. // matrix by diagonal matrix -> matrix operations
  1594.  
  1595. Matrix
  1596. operator + (const Matrix& m, const DiagMatrix& a)
  1597. {
  1598.   int nr = m.rows ();
  1599.   int nc = m.cols ();
  1600.  
  1601.   int a_nr = a.rows ();
  1602.   int a_nc = a.cols ();
  1603.  
  1604.   if (nr != a_nr || nc != a_nc)
  1605.     {
  1606.       gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
  1607.       return Matrix ();
  1608.     }
  1609.  
  1610.   if (nr == 0 || nc == 0)
  1611.     return Matrix (nr, nc);
  1612.  
  1613.   Matrix result (m);
  1614.   int a_len = a.length ();
  1615.   for (int i = 0; i < a_len; i++)
  1616.     result.elem (i, i) += a.elem (i, i);
  1617.  
  1618.   return result;
  1619. }
  1620.  
  1621. Matrix
  1622. operator - (const Matrix& m, const DiagMatrix& a)
  1623. {
  1624.   int nr = m.rows ();
  1625.   int nc = m.cols ();
  1626.  
  1627.   int a_nr = a.rows ();
  1628.   int a_nc = a.cols ();
  1629.  
  1630.   if (nr != a_nr || nc != a_nc)
  1631.     {
  1632.       gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
  1633.       return Matrix ();
  1634.     }
  1635.  
  1636.   if (nr == 0 || nc == 0)
  1637.     return Matrix (nr, nc);
  1638.  
  1639.   Matrix result (m);
  1640.   int a_len = a.length ();
  1641.   for (int i = 0; i < a_len; i++)
  1642.     result.elem (i, i) -= a.elem (i, i);
  1643.  
  1644.   return result;
  1645. }
  1646.  
  1647. Matrix
  1648. operator * (const Matrix& m, const DiagMatrix& a)
  1649. {
  1650.   Matrix retval;
  1651.  
  1652.   int nr = m.rows ();
  1653.   int nc = m.cols ();
  1654.  
  1655.   int a_nr = a.rows ();
  1656.   int a_nc = a.cols ();
  1657.  
  1658.   if (nc != a_nr)
  1659.     gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  1660.   else
  1661.     {
  1662.       if (nr == 0 || nc == 0 || a_nc == 0)
  1663.     retval.resize (nr, a_nc, 0.0);
  1664.       else
  1665.     {
  1666.       retval.resize (nr, a_nc);
  1667.       double *c = retval.fortran_vec ();
  1668.  
  1669.       double *ctmp = 0;
  1670.  
  1671.       int a_len = a.length ();
  1672.  
  1673.       for (int j = 0; j < a_len; j++)
  1674.         {
  1675.           int idx = j * nr;
  1676.           ctmp = c + idx;
  1677.  
  1678.           if (a.elem (j, j) == 1.0)
  1679.         {
  1680.           for (int i = 0; i < nr; i++)
  1681.             ctmp[i] = m.elem (i, j);
  1682.         }
  1683.           else if (a.elem (j, j) == 0.0)
  1684.         {
  1685.           for (int i = 0; i < nr; i++)
  1686.             ctmp[i] = 0.0;
  1687.         }
  1688.           else
  1689.         {
  1690.           for (int i = 0; i < nr; i++)
  1691.             ctmp[i] = a.elem (j, j) * m.elem (i, j);
  1692.         }
  1693.         }
  1694.  
  1695.       if (a_nr < a_nc)
  1696.         {
  1697.           for (int i = nr * nc; i < nr * a_nc; i++)
  1698.         ctmp[i] = 0.0;
  1699.         }
  1700.     }
  1701.     }
  1702.  
  1703.   return retval;
  1704. }
  1705.  
  1706. // diagonal matrix by matrix -> matrix operations
  1707.  
  1708. Matrix
  1709. operator + (const DiagMatrix& m, const Matrix& a)
  1710. {
  1711.   int nr = m.rows ();
  1712.   int nc = m.cols ();
  1713.  
  1714.   int a_nr = a.rows ();
  1715.   int a_nc = a.cols ();
  1716.  
  1717.   if (nr != a_nr || nc != a_nc)
  1718.     {
  1719.       gripe_nonconformant ("operator +", nr, nc, a_nr, a_nc);
  1720.       return Matrix ();
  1721.     }
  1722.  
  1723.   if (nr == 0 || nc == 0)
  1724.     return Matrix (nr, nc);
  1725.  
  1726.   Matrix result (a);
  1727.   for (int i = 0; i < m.length (); i++)
  1728.     result.elem (i, i) += m.elem (i, i);
  1729.  
  1730.   return result;
  1731. }
  1732.  
  1733. Matrix
  1734. operator - (const DiagMatrix& m, const Matrix& a)
  1735. {
  1736.   int nr = m.rows ();
  1737.   int nc = m.cols ();
  1738.  
  1739.   int a_nr = a.rows ();
  1740.   int a_nc = a.cols ();
  1741.  
  1742.   if (nr != a_nr || nc != a_nc)
  1743.     {
  1744.       gripe_nonconformant ("operator -", nr, nc, a_nr, a_nc);
  1745.       return Matrix ();
  1746.     }
  1747.  
  1748.   if (nr == 0 || nc == 0)
  1749.     return Matrix (nr, nc);
  1750.  
  1751.   Matrix result (-a);
  1752.   for (int i = 0; i < m.length (); i++)
  1753.     result.elem (i, i) += m.elem (i, i);
  1754.  
  1755.   return result;
  1756. }
  1757.  
  1758. Matrix
  1759. operator * (const DiagMatrix& m, const Matrix& a)
  1760. {
  1761.   int nr = m.rows ();
  1762.   int nc = m.cols ();
  1763.   int a_nr = a.rows ();
  1764.   int a_nc = a.cols ();
  1765.   if (nc != a_nr)
  1766.     {
  1767.       gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  1768.       return Matrix ();
  1769.     }
  1770.  
  1771.   if (nr == 0 || nc == 0 || a_nc == 0)
  1772.     return Matrix (nr, a_nc, 0.0);
  1773.  
  1774.   Matrix c (nr, a_nc);
  1775.  
  1776.   for (int i = 0; i < m.length (); i++)
  1777.     {
  1778.       if (m.elem (i, i) == 1.0)
  1779.     {
  1780.       for (int j = 0; j < a_nc; j++)
  1781.         c.elem (i, j) = a.elem (i, j);
  1782.     }
  1783.       else if (m.elem (i, i) == 0.0)
  1784.     {
  1785.       for (int j = 0; j < a_nc; j++)
  1786.         c.elem (i, j) = 0.0;
  1787.     }
  1788.       else
  1789.     {
  1790.       for (int j = 0; j < a_nc; j++)
  1791.         c.elem (i, j) = m.elem (i, i) * a.elem (i, j);
  1792.     }
  1793.     }
  1794.  
  1795.   if (nr > nc)
  1796.     {
  1797.       for (int j = 0; j < a_nc; j++)
  1798.     for (int i = a_nr; i < nr; i++)
  1799.       c.elem (i, j) = 0.0;
  1800.     }
  1801.  
  1802.   return c;
  1803. }
  1804.  
  1805. // matrix by matrix -> matrix operations
  1806.  
  1807. Matrix
  1808. operator * (const Matrix& m, const Matrix& a)
  1809. {
  1810.   Matrix retval;
  1811.  
  1812.   int nr = m.rows ();
  1813.   int nc = m.cols ();
  1814.  
  1815.   int a_nr = a.rows ();
  1816.   int a_nc = a.cols ();
  1817.  
  1818.   if (nc != a_nr)
  1819.     gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
  1820.   else
  1821.     {
  1822.       if (nr == 0 || nc == 0 || a_nc == 0)
  1823.     retval.resize (nr, a_nc, 0.0);
  1824.       else
  1825.     {
  1826.       int ld  = nr;
  1827.       int lda = a_nr;
  1828.  
  1829.       retval.resize (nr, a_nc);
  1830.       double *c = retval.fortran_vec ();
  1831.  
  1832.       F77_XFCN (dgemm, DGEMM, ("N", "N", nr, a_nc, nc, 1.0,
  1833.                    m.data (), ld, a.data (), lda, 0.0,
  1834.                    c, nr, 1L, 1L));
  1835.  
  1836.       if (f77_exception_encountered)
  1837.         (*current_liboctave_error_handler)
  1838.           ("unrecoverable error in dgemm");
  1839.     }
  1840.     }
  1841.  
  1842.   return retval;
  1843. }
  1844.  
  1845. // other operations.
  1846.  
  1847. Matrix
  1848. Matrix::map (d_d_Mapper f) const
  1849. {
  1850.   Matrix b (*this);
  1851.   return b.apply (f);
  1852. }
  1853.  
  1854. Matrix&
  1855. Matrix::apply (d_d_Mapper f)
  1856. {
  1857.   double *d = fortran_vec (); // Ensures only one reference to my privates!
  1858.  
  1859.   for (int i = 0; i < length (); i++)
  1860.     d[i] = f (d[i]);
  1861.  
  1862.   return *this;
  1863. }
  1864.  
  1865. bool
  1866. Matrix::any_element_is_negative (void) const
  1867. {
  1868.   int nr = rows ();
  1869.   int nc = cols ();
  1870.  
  1871.   for (int j = 0; j < nc; j++)
  1872.     for (int i = 0; i < nr; i++)
  1873.       if (elem (i, j) < 0.0)
  1874.     return true;
  1875.  
  1876.   return false;
  1877. }
  1878.  
  1879.  
  1880. bool
  1881. Matrix::any_element_is_inf_or_nan (void) const
  1882. {
  1883.   int nr = rows ();
  1884.   int nc = cols ();
  1885.  
  1886.   for (int j = 0; j < nc; j++)
  1887.     for (int i = 0; i < nr; i++)
  1888.       {
  1889.     double val = elem (i, j);
  1890.     if (xisinf (val) || xisnan (val))
  1891.       return 1;
  1892.       }
  1893.   return 0;
  1894. }
  1895.  
  1896. bool
  1897. Matrix::all_elements_are_int_or_inf_or_nan (void) const
  1898. {
  1899.   int nr = rows ();
  1900.   int nc = cols ();
  1901.  
  1902.   for (int j = 0; j < nc; j++)
  1903.     for (int i = 0; i < nr; i++)
  1904.       {
  1905.     double val = elem (i, j);
  1906.     if (xisnan (val) || D_NINT (val) == val)
  1907.       continue;
  1908.     else
  1909.       return false;
  1910.       }
  1911.  
  1912.   return true;
  1913. }
  1914.  
  1915. // Return nonzero if any element of M is not an integer.  Also extract
  1916. // the largest and smallest values and return them in MAX_VAL and MIN_VAL.
  1917.  
  1918. bool
  1919. Matrix::all_integers (double& max_val, double& min_val) const
  1920. {
  1921.   int nr = rows ();
  1922.   int nc = cols ();
  1923.  
  1924.   if (nr > 0 && nc > 0)
  1925.     {
  1926.       max_val = elem (0, 0);
  1927.       min_val = elem (0, 0);
  1928.     }
  1929.   else
  1930.     return false;
  1931.  
  1932.   for (int j = 0; j < nc; j++)
  1933.     for (int i = 0; i < nr; i++)
  1934.       {
  1935.     double val = elem (i, j);
  1936.  
  1937.     if (val > max_val)
  1938.       max_val = val;
  1939.  
  1940.     if (val < min_val)
  1941.       min_val = val;
  1942.  
  1943.     if (D_NINT (val) != val)
  1944.       return false;
  1945.       }
  1946.  
  1947.   return true;
  1948. }
  1949.  
  1950. bool
  1951. Matrix::too_large_for_float (void) const
  1952. {
  1953.   int nr = rows ();
  1954.   int nc = cols ();
  1955.  
  1956.   for (int j = 0; j < nc; j++)
  1957.     for (int i = 0; i < nr; i++)
  1958.       {
  1959.     double val = elem (i, j);
  1960.  
  1961.     if (val > FLT_MAX || val < FLT_MIN)
  1962.       return true;
  1963.       }
  1964.  
  1965.   return false;
  1966. }
  1967.  
  1968. // XXX FIXME XXX Do these really belong here?  They should maybe be
  1969. // cleaned up a bit, no?  What about corresponding functions for the
  1970. // Vectors?
  1971.  
  1972. Matrix
  1973. Matrix::all (void) const
  1974. {
  1975.   int nr = rows ();
  1976.   int nc = cols ();
  1977.   Matrix retval;
  1978.   if (nr > 0 && nc > 0)
  1979.     {
  1980.       if (nr == 1)
  1981.     {
  1982.       retval.resize (1, 1);
  1983.       retval.elem (0, 0) = 1.0;
  1984.       for (int j = 0; j < nc; j++)
  1985.         {
  1986.           if (elem (0, j) == 0.0)
  1987.         {
  1988.           retval.elem (0, 0) = 0.0;
  1989.           break;
  1990.         }
  1991.         }
  1992.     }
  1993.       else if (nc == 1)
  1994.     {
  1995.       retval.resize (1, 1);
  1996.       retval.elem (0, 0) = 1.0;
  1997.       for (int i = 0; i < nr; i++)
  1998.         {
  1999.           if (elem (i, 0) == 0.0)
  2000.         {
  2001.           retval.elem (0, 0) = 0.0;
  2002.           break;
  2003.         }
  2004.         }
  2005.     }
  2006.       else
  2007.     {
  2008.       retval.resize (1, nc);
  2009.       for (int j = 0; j < nc; j++)
  2010.         {
  2011.           retval.elem (0, j) = 1.0;
  2012.           for (int i = 0; i < nr; i++)
  2013.         {
  2014.           if (elem (i, j) == 0.0)
  2015.             {
  2016.               retval.elem (0, j) = 0.0;
  2017.               break;
  2018.             }
  2019.         }
  2020.         }
  2021.     }
  2022.     }
  2023.   return retval;
  2024. }
  2025.  
  2026. Matrix
  2027. Matrix::any (void) const
  2028. {
  2029.   int nr = rows ();
  2030.   int nc = cols ();
  2031.   Matrix retval;
  2032.   if (nr > 0 && nc > 0)
  2033.     {
  2034.       if (nr == 1)
  2035.     {
  2036.       retval.resize (1, 1);
  2037.       retval.elem (0, 0) = 0.0;
  2038.       for (int j = 0; j < nc; j++)
  2039.         {
  2040.           if (elem (0, j) != 0.0)
  2041.         {
  2042.           retval.elem (0, 0) = 1.0;
  2043.           break;
  2044.         }
  2045.         }
  2046.     }
  2047.       else if (nc == 1)
  2048.     {
  2049.       retval.resize (1, 1);
  2050.       retval.elem (0, 0) = 0.0;
  2051.       for (int i = 0; i < nr; i++)
  2052.         {
  2053.           if (elem (i, 0) != 0.0)
  2054.         {
  2055.           retval.elem (0, 0) = 1.0;
  2056.           break;
  2057.         }
  2058.         }
  2059.     }
  2060.       else
  2061.     {
  2062.       retval.resize (1, nc);
  2063.       for (int j = 0; j < nc; j++)
  2064.         {
  2065.           retval.elem (0, j) = 0.0;
  2066.           for (int i = 0; i < nr; i++)
  2067.         {
  2068.           if (elem (i, j) != 0.0)
  2069.             {
  2070.               retval.elem (0, j) = 1.0;
  2071.               break;
  2072.             }
  2073.         }
  2074.         }
  2075.     }
  2076.     }
  2077.   return retval;
  2078. }
  2079.  
  2080. Matrix
  2081. Matrix::cumprod (void) const
  2082. {
  2083.   Matrix retval;
  2084.  
  2085.   int nr = rows ();
  2086.   int nc = cols ();
  2087.  
  2088.   if (nr == 1)
  2089.     {
  2090.       retval.resize (1, nc);
  2091.       if (nc > 0)
  2092.     {
  2093.       double prod = elem (0, 0);
  2094.       for (int j = 0; j < nc; j++)
  2095.         {
  2096.           retval.elem (0, j) = prod;
  2097.           if (j < nc - 1)
  2098.         prod *= elem (0, j+1);
  2099.         }
  2100.     }
  2101.     }
  2102.   else if (nc == 1)
  2103.     {
  2104.       retval.resize (nr, 1);
  2105.       if (nr > 0)
  2106.     {
  2107.       double prod = elem (0, 0);
  2108.       for (int i = 0; i < nr; i++)
  2109.         {
  2110.           retval.elem (i, 0) = prod;
  2111.           if (i < nr - 1)
  2112.         prod *= elem (i+1, 0);
  2113.         }
  2114.     }
  2115.     }
  2116.   else
  2117.     {
  2118.       retval.resize (nr, nc);
  2119.       if (nr > 0 && nc > 0)
  2120.     {
  2121.       for (int j = 0; j < nc; j++)
  2122.         {
  2123.           double prod = elem (0, j);
  2124.           for (int i = 0; i < nr; i++)
  2125.         {
  2126.           retval.elem (i, j) = prod;
  2127.           if (i < nr - 1)
  2128.             prod *= elem (i+1, j);
  2129.         }
  2130.         }
  2131.     }
  2132.     }
  2133.   return retval;
  2134. }
  2135.  
  2136. Matrix
  2137. Matrix::cumsum (void) const
  2138. {
  2139.   Matrix retval;
  2140.  
  2141.   int nr = rows ();
  2142.   int nc = cols ();
  2143.  
  2144.   if (nr == 1)
  2145.     {
  2146.       retval.resize (1, nc);
  2147.       if (nc > 0)
  2148.     {
  2149.       double sum = elem (0, 0);
  2150.       for (int j = 0; j < nc; j++)
  2151.         {
  2152.           retval.elem (0, j) = sum;
  2153.           if (j < nc - 1)
  2154.         sum += elem (0, j+1);
  2155.         }
  2156.     }
  2157.     }
  2158.   else if (nc == 1)
  2159.     {
  2160.       retval.resize (nr, 1);
  2161.       if (nr > 0)
  2162.     {
  2163.       double sum = elem (0, 0);
  2164.       for (int i = 0; i < nr; i++)
  2165.         {
  2166.           retval.elem (i, 0) = sum;
  2167.           if (i < nr - 1)
  2168.         sum += elem (i+1, 0);
  2169.         }
  2170.     }
  2171.     }
  2172.   else
  2173.     {
  2174.       retval.resize (nr, nc);
  2175.       if (nr > 0 && nc > 0)
  2176.     {
  2177.       for (int j = 0; j < nc; j++)
  2178.         {
  2179.           double sum = elem (0, j);
  2180.           for (int i = 0; i < nr; i++)
  2181.         {
  2182.           retval.elem (i, j) = sum;
  2183.           if (i < nr - 1)
  2184.             sum += elem (i+1, j);
  2185.         }
  2186.         }
  2187.     }
  2188.     }
  2189.   return retval;
  2190. }
  2191.  
  2192. Matrix
  2193. Matrix::prod (void) const
  2194. {
  2195.   Matrix retval;
  2196.  
  2197.   int nr = rows ();
  2198.   int nc = cols ();
  2199.  
  2200.   if (nr == 1)
  2201.     {
  2202.       retval.resize (1, 1);
  2203.       retval.elem (0, 0) = 1.0;
  2204.       for (int j = 0; j < nc; j++)
  2205.     retval.elem (0, 0) *= elem (0, j);
  2206.     }
  2207.   else if (nc == 1)
  2208.     {
  2209.       retval.resize (1, 1);
  2210.       retval.elem (0, 0) = 1.0;
  2211.       for (int i = 0; i < nr; i++)
  2212.     retval.elem (0, 0) *= elem (i, 0);
  2213.     }
  2214.   else
  2215.     {
  2216.       if (nc == 0)
  2217.     {
  2218.       retval.resize (1, 1);
  2219.       retval.elem (0, 0) = 1.0;
  2220.     }
  2221.       else
  2222.     retval.resize (1, nc);
  2223.  
  2224.       for (int j = 0; j < nc; j++)
  2225.     {
  2226.       retval.elem (0, j) = 1.0;
  2227.       for (int i = 0; i < nr; i++)
  2228.         retval.elem (0, j) *= elem (i, j);
  2229.     }
  2230.     }
  2231.   return retval;
  2232. }
  2233.  
  2234. Matrix
  2235. Matrix::sum (void) const
  2236. {
  2237.   Matrix retval;
  2238.  
  2239.   int nr = rows ();
  2240.   int nc = cols ();
  2241.  
  2242.   if (nr == 1)
  2243.     {
  2244.       retval.resize (1, 1);
  2245.       retval.elem (0, 0) = 0.0;
  2246.       for (int j = 0; j < nc; j++)
  2247.     retval.elem (0, 0) += elem (0, j);
  2248.     }
  2249.   else if (nc == 1)
  2250.     {
  2251.       retval.resize (1, 1);
  2252.       retval.elem (0, 0) = 0.0;
  2253.       for (int i = 0; i < nr; i++)
  2254.     retval.elem (0, 0) += elem (i, 0);
  2255.     }
  2256.   else
  2257.     {
  2258.       if (nc == 0)
  2259.     {
  2260.       retval.resize (1, 1);
  2261.       retval.elem (0, 0) = 0.0;
  2262.     }
  2263.       else
  2264.     retval.resize (1, nc);
  2265.  
  2266.       for (int j = 0; j < nc; j++)
  2267.     {
  2268.       retval.elem (0, j) = 0.0;
  2269.       for (int i = 0; i < nr; i++)
  2270.         retval.elem (0, j) += elem (i, j);
  2271.     }
  2272.     }
  2273.   return retval;
  2274. }
  2275.  
  2276. Matrix
  2277. Matrix::sumsq (void) const
  2278. {
  2279.   Matrix retval;
  2280.  
  2281.   int nr = rows ();
  2282.   int nc = cols ();
  2283.  
  2284.   if (nr == 1)
  2285.     {
  2286.       retval.resize (1, 1);
  2287.       retval.elem (0, 0) = 0.0;
  2288.       for (int j = 0; j < nc; j++)
  2289.     {
  2290.       double d = elem (0, j);
  2291.       retval.elem (0, 0) += d * d;
  2292.     }
  2293.     }
  2294.   else if (nc == 1)
  2295.     {
  2296.       retval.resize (1, 1);
  2297.       retval.elem (0, 0) = 0.0;
  2298.       for (int i = 0; i < nr; i++)
  2299.     {
  2300.       double d = elem (i, 0);
  2301.       retval.elem (0, 0) += d * d;
  2302.     }
  2303.     }
  2304.   else
  2305.     {
  2306.       retval.resize (1, nc);
  2307.       for (int j = 0; j < nc; j++)
  2308.     {
  2309.       retval.elem (0, j) = 0.0;
  2310.       for (int i = 0; i < nr; i++)
  2311.         {
  2312.           double d = elem (i, j);
  2313.           retval.elem (0, j) += d * d;
  2314.         }
  2315.     }
  2316.     }
  2317.   return retval;
  2318. }
  2319.  
  2320. Matrix
  2321. Matrix::abs (void) const
  2322. {
  2323.   int nr = rows ();
  2324.   int nc = cols ();
  2325.  
  2326.   Matrix retval (nr, nc);
  2327.  
  2328.   for (int j = 0; j < nc; j++)
  2329.     for (int i = 0; i < nr; i++)
  2330.       retval (i, j) = fabs (elem (i, j));
  2331.  
  2332.   return retval;
  2333. }
  2334.  
  2335. ColumnVector
  2336. Matrix::diag (void) const
  2337. {
  2338.   return diag (0);
  2339. }
  2340.  
  2341. ColumnVector
  2342. Matrix::diag (int k) const
  2343. {
  2344.   int nnr = rows ();
  2345.   int nnc = cols ();
  2346.   if (k > 0)
  2347.     nnc -= k;
  2348.   else if (k < 0)
  2349.     nnr += k;
  2350.  
  2351.   ColumnVector d;
  2352.  
  2353.   if (nnr > 0 && nnc > 0)
  2354.     {
  2355.       int ndiag = (nnr < nnc) ? nnr : nnc;
  2356.  
  2357.       d.resize (ndiag);
  2358.  
  2359.       if (k > 0)
  2360.     {
  2361.       for (int i = 0; i < ndiag; i++)
  2362.         d.elem (i) = elem (i, i+k);
  2363.     }
  2364.       else if ( k < 0)
  2365.     {
  2366.       for (int i = 0; i < ndiag; i++)
  2367.         d.elem (i) = elem (i-k, i);
  2368.     }
  2369.       else
  2370.     {
  2371.       for (int i = 0; i < ndiag; i++)
  2372.         d.elem (i) = elem (i, i);
  2373.     }
  2374.     }
  2375.   else
  2376.     cerr << "diag: requested diagonal out of range\n";
  2377.  
  2378.   return d;
  2379. }
  2380.  
  2381. ColumnVector
  2382. Matrix::row_min (void) const
  2383. {
  2384.   Array<int> index;
  2385.   return row_min (index);
  2386. }
  2387.  
  2388. ColumnVector
  2389. Matrix::row_min (Array<int>& index) const
  2390. {
  2391.   ColumnVector result;
  2392.  
  2393.   int nr = rows ();
  2394.   int nc = cols ();
  2395.  
  2396.   if (nr > 0 && nc > 0)
  2397.     {
  2398.       result.resize (nr);
  2399.       index.resize (nr);
  2400.  
  2401.       for (int i = 0; i < nr; i++)
  2402.         {
  2403.       int idx = 0;
  2404.  
  2405.       double tmp_min = elem (i, idx);
  2406.  
  2407.       if (xisnan (tmp_min))
  2408.         idx = -1;
  2409.       else
  2410.         {
  2411.           for (int j = 1; j < nc; j++)
  2412.         {
  2413.           double tmp = elem (i, j);
  2414.  
  2415.           if (xisnan (tmp))
  2416.             {
  2417.               idx = -1;
  2418.               break;
  2419.             }
  2420.           else if (tmp < tmp_min)
  2421.             {
  2422.               idx = j;
  2423.               tmp_min = tmp;
  2424.             }
  2425.         }
  2426.         }
  2427.  
  2428.       result.elem (i) = (idx < 0) ? octave_NaN : tmp_min;
  2429.       index.elem (i) = idx;
  2430.         }
  2431.     }
  2432.  
  2433.   return result;
  2434. }
  2435.  
  2436. ColumnVector
  2437. Matrix::row_max (void) const
  2438. {
  2439.   Array<int> index;
  2440.   return row_max (index);
  2441. }
  2442.  
  2443. ColumnVector
  2444. Matrix::row_max (Array<int>& index) const
  2445. {
  2446.   ColumnVector result;
  2447.  
  2448.   int nr = rows ();
  2449.   int nc = cols ();
  2450.  
  2451.   if (nr > 0 && nc > 0)
  2452.     {
  2453.       result.resize (nr);
  2454.       index.resize (nr);
  2455.  
  2456.       for (int i = 0; i < nr; i++)
  2457.         {
  2458.       int idx = 0;
  2459.  
  2460.       double tmp_max = elem (i, idx);
  2461.  
  2462.       if (xisnan (tmp_max))
  2463.         idx = -1;
  2464.       else
  2465.         {
  2466.           for (int j = 1; j < nc; j++)
  2467.         {
  2468.           double tmp = elem (i, j);
  2469.  
  2470.           if (xisnan (tmp))
  2471.             {
  2472.               idx = -1;
  2473.               break;
  2474.             }
  2475.           else if (tmp > tmp_max)
  2476.             {
  2477.               idx = j;
  2478.               tmp_max = tmp;
  2479.             }
  2480.         }
  2481.         }
  2482.  
  2483.       result.elem (i) = (idx < 0) ? octave_NaN : tmp_max;
  2484.       index.elem (i) = idx;
  2485.         }
  2486.     }
  2487.  
  2488.   return result;
  2489. }
  2490.  
  2491. RowVector
  2492. Matrix::column_min (void) const
  2493. {
  2494.   Array<int> index;
  2495.   return column_min (index);
  2496. }
  2497.  
  2498. RowVector
  2499. Matrix::column_min (Array<int>& index) const
  2500. {
  2501.   RowVector result;
  2502.  
  2503.   int nr = rows ();
  2504.   int nc = cols ();
  2505.  
  2506.   if (nr > 0 && nc > 0)
  2507.     {
  2508.       result.resize (nc);
  2509.       index.resize (nc);
  2510.  
  2511.       for (int j = 0; j < nc; j++)
  2512.         {
  2513.       int idx = 0;
  2514.  
  2515.       double tmp_min = elem (idx, j);
  2516.  
  2517.       if (xisnan (tmp_min))
  2518.         idx = -1;
  2519.       else
  2520.         {
  2521.           for (int i = 1; i < nr; i++)
  2522.         {
  2523.           double tmp = elem (i, j);
  2524.  
  2525.           if (xisnan (tmp))
  2526.             {
  2527.               idx = -1;
  2528.               break;
  2529.             }
  2530.           else if (tmp < tmp_min)
  2531.             {
  2532.               idx = i;
  2533.               tmp_min = tmp;
  2534.             }
  2535.         }
  2536.         }
  2537.  
  2538.       result.elem (j) = (idx < 0) ? octave_NaN : tmp_min;
  2539.       index.elem (j) = idx;
  2540.         }
  2541.     }
  2542.  
  2543.   return result;
  2544. }
  2545.  
  2546. RowVector
  2547. Matrix::column_max (void) const
  2548. {
  2549.   Array<int> index;
  2550.   return column_max (index);
  2551. }
  2552.  
  2553. RowVector
  2554. Matrix::column_max (Array<int>& index) const
  2555. {
  2556.   RowVector result;
  2557.  
  2558.   int nr = rows ();
  2559.   int nc = cols ();
  2560.  
  2561.   if (nr > 0 && nc > 0)
  2562.     {
  2563.       result.resize (nc);
  2564.       index.resize (nc);
  2565.  
  2566.       for (int j = 0; j < nc; j++)
  2567.         {
  2568.       int idx = 0;
  2569.  
  2570.       double tmp_max = elem (idx, j);
  2571.  
  2572.       if (xisnan (tmp_max))
  2573.         idx = -1;
  2574.       else
  2575.         {
  2576.           for (int i = 1; i < nr; i++)
  2577.         {
  2578.           double tmp = elem (i, j);
  2579.  
  2580.           if (xisnan (tmp))
  2581.             {
  2582.               idx = -1;
  2583.               break;
  2584.             }
  2585.           else if (tmp > tmp_max)
  2586.             {
  2587.               idx = i;
  2588.               tmp_max = tmp;
  2589.             }
  2590.         }
  2591.         }
  2592.  
  2593.       result.elem (j) = (idx < 0) ? octave_NaN : tmp_max;
  2594.       index.elem (j) = idx;
  2595.         }
  2596.     }
  2597.  
  2598.   return result;
  2599. }
  2600.  
  2601. ostream&
  2602. operator << (ostream& os, const Matrix& a)
  2603. {
  2604. //  int field_width = os.precision () + 7;
  2605.  
  2606.   for (int i = 0; i < a.rows (); i++)
  2607.     {
  2608.       for (int j = 0; j < a.cols (); j++)
  2609.     os << " " /* setw (field_width) */ << a.elem (i, j);
  2610.       os << "\n";
  2611.     }
  2612.   return os;
  2613. }
  2614.  
  2615. istream&
  2616. operator >> (istream& is, Matrix& a)
  2617. {
  2618.   int nr = a.rows ();
  2619.   int nc = a.cols ();
  2620.  
  2621.   if (nr < 1 || nc < 1)
  2622.     is.clear (ios::badbit);
  2623.   else
  2624.     {
  2625.       double tmp;
  2626.       for (int i = 0; i < nr; i++)
  2627.     for (int j = 0; j < nc; j++)
  2628.       {
  2629.         is >> tmp;
  2630.         if (is)
  2631.           a.elem (i, j) = tmp;
  2632.         else
  2633.           goto done;
  2634.       }
  2635.     }
  2636.  
  2637. done:
  2638.  
  2639.   return is;
  2640. }
  2641.  
  2642. template <class T>
  2643. static void
  2644. read_int (istream& is, bool swap_bytes, T& val)
  2645. {
  2646.   is.read ((char *) &val, sizeof (T));
  2647.  
  2648.   if (swap_bytes)
  2649.     {
  2650.       switch (sizeof (T))
  2651.     {
  2652.     case 1:
  2653.       break;
  2654.  
  2655.     case 2:
  2656.       swap_2_bytes ((char *) &val);
  2657.       break;
  2658.  
  2659.     case 4:
  2660.       swap_4_bytes ((char *) &val);
  2661.       break;
  2662.  
  2663.     case 8:
  2664.       swap_8_bytes ((char *) &val);
  2665.       break;
  2666.  
  2667.     default:
  2668.       (*current_liboctave_error_handler)
  2669.         ("read_int: unrecognized data format!");
  2670.     }
  2671.     }
  2672. }
  2673.  
  2674. template void read_int (istream&, bool, char&);
  2675. template void read_int (istream&, bool, signed char&);
  2676. template void read_int (istream&, bool, unsigned char&);
  2677. template void read_int (istream&, bool, short&);
  2678. template void read_int (istream&, bool, unsigned short&);
  2679. template void read_int (istream&, bool, int&);
  2680. template void read_int (istream&, bool, unsigned int&);
  2681. template void read_int (istream&, bool, long&);
  2682. template void read_int (istream&, bool, unsigned long&);
  2683.  
  2684. static inline bool
  2685. do_read (istream& is, oct_data_conv::data_type dt, 
  2686.      oct_mach_info::float_format flt_fmt, bool swap_bytes,
  2687.      bool do_float_conversion, double& val)
  2688. {
  2689.   bool retval = true;
  2690.  
  2691.   switch (dt)
  2692.     {
  2693.     case oct_data_conv::dt_char:
  2694.       {
  2695.     char tmp;
  2696.     read_int (is, swap_bytes, tmp);
  2697.     val = tmp;
  2698.       }
  2699.       break;
  2700.  
  2701.     case oct_data_conv::dt_schar:
  2702.       {
  2703.     signed char tmp;
  2704.     read_int (is, swap_bytes, tmp);
  2705.     val = tmp;
  2706.       }
  2707.       break;
  2708.  
  2709.     case oct_data_conv::dt_uchar:
  2710.       {
  2711.     unsigned char tmp;
  2712.     read_int (is, swap_bytes, tmp);
  2713.     val = tmp;
  2714.       }
  2715.       break;
  2716.  
  2717.     case oct_data_conv::dt_short:
  2718.       {
  2719.     short tmp;
  2720.     read_int (is, swap_bytes, tmp);
  2721.     val = tmp;
  2722.       }
  2723.       break;
  2724.  
  2725.     case oct_data_conv::dt_ushort:
  2726.       {
  2727.     unsigned short tmp;
  2728.     read_int (is, swap_bytes, tmp);
  2729.     val = tmp;
  2730.       }
  2731.       break;
  2732.  
  2733.     case oct_data_conv::dt_int:
  2734.       {
  2735.     int tmp;
  2736.     read_int (is, swap_bytes, tmp);
  2737.     val = tmp;
  2738.       }
  2739.       break;
  2740.  
  2741.     case oct_data_conv::dt_uint:
  2742.       {
  2743.     unsigned int tmp;
  2744.     read_int (is, swap_bytes, tmp);
  2745.     val = tmp;
  2746.       }
  2747.       break;
  2748.  
  2749.     case oct_data_conv::dt_long:
  2750.       {
  2751.     long tmp;
  2752.     read_int (is, swap_bytes, tmp);
  2753.     val = tmp;
  2754.       }
  2755.       break;
  2756.  
  2757.     case oct_data_conv::dt_ulong:
  2758.       {
  2759.     unsigned long tmp;
  2760.     read_int (is, swap_bytes, tmp);
  2761.     val = tmp;
  2762.       }
  2763.       break;
  2764.  
  2765.     case oct_data_conv::dt_float:
  2766.       {
  2767.     float f;
  2768.  
  2769.     is.read ((char *) &f, sizeof (float));
  2770.  
  2771.     if (do_float_conversion)
  2772.       do_float_format_conversion (&f, 1, flt_fmt);
  2773.  
  2774.     val = f;
  2775.       }
  2776.       break;
  2777.  
  2778.     case oct_data_conv::dt_double:
  2779.       {
  2780.     is.read ((char *) &val, sizeof (double));
  2781.  
  2782.     if (do_float_conversion)
  2783.       do_double_format_conversion (&val, 1, flt_fmt);
  2784.       }
  2785.       break;
  2786.  
  2787.     default:
  2788.       retval = false;
  2789.       (*current_liboctave_error_handler)
  2790.     ("read: invalid type specification");
  2791.       break;
  2792.     }
  2793.  
  2794.   return retval;
  2795. }
  2796.  
  2797. int
  2798. Matrix::read (istream& is, int nr, int nc,
  2799.           oct_data_conv::data_type dt, int skip,
  2800.           oct_mach_info::float_format flt_fmt)
  2801. {
  2802.   int retval = -1;
  2803.  
  2804.   bool ok = true;
  2805.  
  2806.   int count = 0;
  2807.  
  2808.   double *data = 0;
  2809.   int max_size = 0;
  2810.  
  2811.   int final_nr = 0;
  2812.   int final_nc = 0;
  2813.  
  2814.   if (nr > 0)
  2815.     {
  2816.       if (nc > 0)
  2817.     {
  2818.       resize (nr, nc, 0.0);
  2819.       data = fortran_vec ();
  2820.       max_size = nr * nc;
  2821.     }
  2822.       else
  2823.     {
  2824.       resize (nr, 32, 0.0);
  2825.       data = fortran_vec ();
  2826.       max_size = nr * 32;
  2827.     }
  2828.     }
  2829.   else
  2830.     {
  2831.       resize (32, 1, 0.0);
  2832.       data = fortran_vec ();
  2833.       max_size = 32;
  2834.     }
  2835.  
  2836.   oct_mach_info::float_format native_flt_fmt
  2837.     = oct_mach_info::float_format ();
  2838.  
  2839.   bool do_float_conversion = (flt_fmt != native_flt_fmt);
  2840.  
  2841.   // XXX FIXME XXX -- byte order for Cray?
  2842.  
  2843.   bool swap_bytes = false;
  2844.  
  2845.   if (oct_mach_info::words_big_endian ())
  2846.     swap_bytes = (flt_fmt == oct_mach_info::ieee_little_endian
  2847.          || flt_fmt == oct_mach_info::vax_g
  2848.          || flt_fmt == oct_mach_info::vax_g);
  2849.   else
  2850.     swap_bytes = (flt_fmt == oct_mach_info::ieee_big_endian);
  2851.  
  2852.   for (;;)
  2853.     {
  2854.       // XXX FIXME XXX -- maybe there should be a special case for
  2855.       // skip == 0.
  2856.  
  2857.       if (is)
  2858.     {
  2859.       if (nr > 0 && nc > 0 && count == max_size)
  2860.         {
  2861.           final_nr = nr;
  2862.           final_nc = nc;
  2863.  
  2864.           break;
  2865.         }
  2866.  
  2867.       if (skip != 0)
  2868.         is.seekg (skip, ios::cur);
  2869.  
  2870.       if (is)
  2871.         {
  2872.           double tmp = 0.0;
  2873.  
  2874.           ok = do_read (is, dt, flt_fmt, swap_bytes,
  2875.                 do_float_conversion, tmp);
  2876.  
  2877.           if (ok)
  2878.         {
  2879.           if (is)
  2880.             {
  2881.               if (count == max_size)
  2882.             {
  2883.               max_size *= 2;
  2884.  
  2885.               if (nr > 0)
  2886.                 resize (nr, max_size / nr, 0.0);
  2887.               else
  2888.                 resize (max_size, 1, 0.0);
  2889.  
  2890.               data = fortran_vec ();
  2891.             }
  2892.  
  2893.               data[count++] = tmp;
  2894.             }
  2895.           else
  2896.             {
  2897.               if (is.eof ())
  2898.             {
  2899.               if (nr > 0)
  2900.                 {
  2901.                   if (count > nr)
  2902.                 {
  2903.                   final_nr = nr;
  2904.                   final_nc = (count - 1) / nr + 1;
  2905.                 }
  2906.                   else
  2907.                 {
  2908.                   final_nr = count;
  2909.                   final_nc = 1;
  2910.                 }
  2911.                 }
  2912.               else
  2913.                 {
  2914.                   final_nr = count;
  2915.                   final_nc = 1;
  2916.                 }
  2917.             }
  2918.  
  2919.               break;
  2920.             }
  2921.         }
  2922.           else
  2923.         break;
  2924.         }
  2925.       else
  2926.         {
  2927.           ok = false;
  2928.           break;
  2929.         }
  2930.     }
  2931.       else
  2932.     {
  2933.       ok = false;
  2934.       break;
  2935.     }
  2936.     }
  2937.  
  2938.   if (ok)
  2939.     {
  2940.       resize (final_nr, final_nc, 0.0);
  2941.  
  2942.       retval = count;
  2943.     }
  2944.  
  2945.   return retval;
  2946. }
  2947.  
  2948. template <class T>
  2949. static void
  2950. write_int (ostream& os, bool swap_bytes, T val)
  2951. {
  2952.   if (swap_bytes)
  2953.     {
  2954.       switch (sizeof (T))
  2955.     {
  2956.     case 1:
  2957.       break;
  2958.  
  2959.     case 2:
  2960.       swap_2_bytes ((char *) &val);
  2961.       break;
  2962.  
  2963.     case 4:
  2964.       swap_4_bytes ((char *) &val);
  2965.       break;
  2966.  
  2967.     case 8:
  2968.       swap_8_bytes ((char *) &val);
  2969.       break;
  2970.  
  2971.     default:
  2972.       (*current_liboctave_error_handler)
  2973.         ("write_int: unrecognized data format!");
  2974.     }
  2975.     }
  2976.  
  2977.   os.write ((char *) &val, sizeof (T));
  2978. }
  2979.  
  2980. template void write_int (ostream&, bool, char);
  2981. template void write_int (ostream&, bool, signed char);
  2982. template void write_int (ostream&, bool, unsigned char);
  2983. template void write_int (ostream&, bool, short);
  2984. template void write_int (ostream&, bool, unsigned short);
  2985. template void write_int (ostream&, bool, int);
  2986. template void write_int (ostream&, bool, unsigned int);
  2987. template void write_int (ostream&, bool, long);
  2988. template void write_int (ostream&, bool, unsigned long);
  2989.  
  2990. static inline bool
  2991. do_write (ostream& os, double d, oct_data_conv::data_type dt,
  2992.       oct_mach_info::float_format flt_fmt, bool swap_bytes,
  2993.       bool do_float_conversion)
  2994. {
  2995.   bool retval = true;
  2996.  
  2997.   switch (dt)
  2998.     {
  2999.     case oct_data_conv::dt_char:
  3000.       write_int (os, swap_bytes, (char) d);
  3001.       break;
  3002.  
  3003.     case oct_data_conv::dt_schar:
  3004.       write_int (os, swap_bytes, (signed char) d);
  3005.       break;
  3006.  
  3007.     case oct_data_conv::dt_uchar:
  3008.       write_int (os, swap_bytes, (unsigned char) d);
  3009.       break;
  3010.  
  3011.     case oct_data_conv::dt_short:
  3012.       write_int (os, swap_bytes, (short) d);
  3013.       break;
  3014.  
  3015.     case oct_data_conv::dt_ushort:
  3016.       write_int (os, swap_bytes, (unsigned short) d);
  3017.       break;
  3018.  
  3019.     case oct_data_conv::dt_int:
  3020.       write_int (os, swap_bytes, (int) d);
  3021.       break;
  3022.  
  3023.     case oct_data_conv::dt_uint:
  3024.       write_int (os, swap_bytes, (unsigned int) d);
  3025.       break;
  3026.  
  3027.     case oct_data_conv::dt_long:
  3028.       write_int (os, swap_bytes, (long) d);
  3029.       break;
  3030.  
  3031.     case oct_data_conv::dt_ulong:
  3032.       write_int (os, swap_bytes, (unsigned long) d);
  3033.       break;
  3034.  
  3035.     case oct_data_conv::dt_float:
  3036.       {
  3037.     float f = (float) d;
  3038.  
  3039.     if (do_float_conversion)
  3040.       do_float_format_conversion (&f, 1, flt_fmt);
  3041.  
  3042.     os.write ((char *) &f, sizeof (float));
  3043.       }
  3044.       break;
  3045.  
  3046.     case oct_data_conv::dt_double:
  3047.       {
  3048.     if (do_float_conversion)
  3049.       do_double_format_conversion (&d, 1, flt_fmt);
  3050.  
  3051.     os.write ((char *) &d, sizeof (double));
  3052.       }
  3053.       break;
  3054.  
  3055.     default:
  3056.       retval = false;
  3057.       (*current_liboctave_error_handler)
  3058.     ("write: invalid type specification");
  3059.       break;
  3060.     }
  3061.  
  3062.   return retval;
  3063. }
  3064.  
  3065. int
  3066. Matrix::write (ostream& os, oct_data_conv::data_type dt, int skip,
  3067.            oct_mach_info::float_format flt_fmt)
  3068. {
  3069.   int retval = -1;
  3070.  
  3071.   bool ok = true;
  3072.  
  3073.   int count = 0;
  3074.  
  3075.   const double *d = data ();
  3076.  
  3077.   int n = length ();
  3078.  
  3079.   oct_mach_info::float_format native_flt_fmt
  3080.     = oct_mach_info::float_format ();
  3081.  
  3082.   bool do_float_conversion = (flt_fmt != native_flt_fmt);
  3083.  
  3084.   // XXX FIXME XXX -- byte order for Cray?
  3085.  
  3086.   bool swap_bytes = false;
  3087.  
  3088.   if (oct_mach_info::words_big_endian ())
  3089.     swap_bytes = (flt_fmt == oct_mach_info::ieee_little_endian
  3090.          || flt_fmt == oct_mach_info::vax_g
  3091.          || flt_fmt == oct_mach_info::vax_g);
  3092.   else
  3093.     swap_bytes = (flt_fmt == oct_mach_info::ieee_big_endian);
  3094.  
  3095.   for (int i = 0; i < n; i++)
  3096.     {
  3097.       if (os)
  3098.     {
  3099.       if (skip != 0)
  3100.         os.seekp (skip, ios::cur);
  3101.  
  3102.       if (os)
  3103.         {
  3104.           ok = do_write (os, d[i], dt, flt_fmt, swap_bytes,
  3105.                  do_float_conversion);
  3106.  
  3107.           if (os && ok)
  3108.         count++;
  3109.           else
  3110.         break;
  3111.         }
  3112.       else
  3113.         {
  3114.           ok = false;
  3115.           break;
  3116.         }
  3117.     }
  3118.       else
  3119.     {
  3120.       ok = false;
  3121.       break;
  3122.     }
  3123.     }
  3124.  
  3125.   if (ok)
  3126.     retval = count;
  3127.  
  3128.   return retval;
  3129. }
  3130.  
  3131.  
  3132.  
  3133. Matrix
  3134. Givens (double x, double y)
  3135. {
  3136.   double cc, s, temp_r;
  3137.  
  3138.   F77_FCN (dlartg, DLARTG) (x, y, cc, s, temp_r);
  3139.  
  3140.   Matrix g (2, 2);
  3141.  
  3142.   g.elem (0, 0) = cc;
  3143.   g.elem (1, 1) = cc;
  3144.   g.elem (0, 1) = s;
  3145.   g.elem (1, 0) = -s;
  3146.  
  3147.   return g;
  3148. }
  3149.  
  3150. Matrix
  3151. Sylvester (const Matrix& a, const Matrix& b, const Matrix& c)
  3152. {
  3153.   Matrix retval;
  3154.  
  3155.   // XXX FIXME XXX -- need to check that a, b, and c are all the same
  3156.   // size.
  3157.  
  3158.   // Compute Schur decompositions.
  3159.  
  3160.   SCHUR as (a, "U");
  3161.   SCHUR bs (b, "U");
  3162.   
  3163.   // Transform c to new coordinates.
  3164.  
  3165.   Matrix ua = as.unitary_matrix ();
  3166.   Matrix sch_a = as.schur_matrix ();
  3167.  
  3168.   Matrix ub = bs.unitary_matrix ();
  3169.   Matrix sch_b = bs.schur_matrix ();
  3170.   
  3171.   Matrix cx = ua.transpose () * c * ub;
  3172.   
  3173.   // Solve the sylvester equation, back-transform, and return the
  3174.   // solution.
  3175.  
  3176.   int a_nr = a.rows ();
  3177.   int b_nr = b.rows ();
  3178.  
  3179.   double scale;
  3180.   int info;
  3181.  
  3182.   double *pa = sch_a.fortran_vec ();
  3183.   double *pb = sch_b.fortran_vec ();
  3184.   double *px = cx.fortran_vec ();
  3185.  
  3186.   F77_XFCN (dtrsyl, DTRSYL, ("N", "N", 1, a_nr, b_nr, pa, a_nr, pb,
  3187.                  b_nr, px, a_nr, scale, info, 1L, 1L));
  3188.  
  3189.  
  3190.   if (f77_exception_encountered)
  3191.     (*current_liboctave_error_handler) ("unrecoverable error in dtrsyl");
  3192.   else
  3193.     {
  3194.       // XXX FIXME XXX -- check info?
  3195.   
  3196.       retval = -ua*cx*ub.transpose ();
  3197.     }
  3198.  
  3199.   return retval;
  3200. }
  3201.  
  3202. ComplexColumnVector
  3203. Qzval (const Matrix& a, const Matrix& b)
  3204. {
  3205.   ComplexColumnVector retval;
  3206.  
  3207.   int a_nr = a.rows();
  3208.   int a_nc = a.cols();
  3209.  
  3210.   int b_nr = b.rows();
  3211.   int b_nc = b.cols();
  3212.  
  3213.   if (a_nr == a_nc)
  3214.     {
  3215.       if (a_nr == b_nr && a_nc == b_nc)
  3216.     {
  3217.       if (a_nr != 0)
  3218.         {
  3219.           Matrix jnk (a_nr, a_nr, 0.0);
  3220.           double *pjnk = jnk.fortran_vec ();
  3221.  
  3222.           ColumnVector alfr (a_nr);
  3223.           double *palfr = alfr.fortran_vec ();
  3224.  
  3225.           ColumnVector alfi (a_nr);
  3226.           double *palfi = alfi.fortran_vec ();
  3227.  
  3228.           ColumnVector beta (a_nr);
  3229.           double *pbeta = beta.fortran_vec ();
  3230.  
  3231.           Matrix atmp = a;
  3232.           double *pa = atmp.fortran_vec ();
  3233.  
  3234.           Matrix btmp = b;
  3235.           double *pb = btmp.fortran_vec ();
  3236.  
  3237.           long matz = 0;
  3238.           int info;
  3239.  
  3240.           // XXX FIXME ??? XXX
  3241.           double eps = DBL_EPSILON;
  3242.  
  3243.           F77_FCN (qzhes, QZHES) (a_nr, a_nr, pa, pb, matz, pjnk);
  3244.  
  3245.           F77_FCN (qzit, QZIT) (a_nr, a_nr, pa, pb, eps, matz, pjnk, info);
  3246.  
  3247.           if (! info)
  3248.         {
  3249.           F77_FCN (qzval, QZVAL) (a_nr, a_nr, pa, pb, palfr,
  3250.                       palfi, pbeta, matz, pjnk);
  3251.  
  3252.           // Count and extract finite generalized eigenvalues.
  3253.  
  3254.           int cnt = 0;
  3255.  
  3256.           for (int i = 0; i < a_nr; i++)
  3257.             if (beta(i) != 0)
  3258.               cnt++;
  3259.  
  3260.           ComplexColumnVector cx (cnt);
  3261.  
  3262.           cnt = 0;
  3263.  
  3264.           for (int i = 0; i < a_nr; i++)
  3265.             {
  3266.               if (beta(i) != 0)
  3267.             {
  3268.               // Finite generalized eigenvalue.
  3269.  
  3270.               cx(cnt++) = Complex (alfr(i), alfi(i)) / beta(i);
  3271.             }
  3272.             }
  3273.  
  3274.           retval = cx;
  3275.         }
  3276.           else
  3277.         (*current_liboctave_error_handler)
  3278.           ("qzval: trouble in qzit, info = %d", info);
  3279.         }
  3280.     }
  3281.       else
  3282.     gripe_nonconformant ("qzval", a_nr, a_nc, b_nr, b_nc);
  3283.     }
  3284.   else
  3285.     (*current_liboctave_error_handler) ("qzval: square matrices required");
  3286.  
  3287.   return retval;
  3288. }
  3289.  
  3290. /*
  3291. ;;; Local Variables: ***
  3292. ;;; mode: C++ ***
  3293. ;;; End: ***
  3294. */
  3295.